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1. Motivation for error estimates of molecular dynamics 

Molecular dynamics is a computational method to study molecular systems in materials science, chemistry 
and molecular biology. The simulations are used, for example, in designing and understanding new materials 
or for determining biochemical reactions in drug design, |14) . The wide popularity of molecular dynamics 
simulations relies on the fact that in many cases it agrees very well with experiments. Indeed when we have 
experimental data it is easy to verify correctness of the method by comparing with experiments at certain 
parameter regimes. However, if we want the simulation to predict something that has no comparing experi- 
ment, we need a mathematical estimate of the accuracy of the computation. In the case of molecular systems 
with few particles such studies are made by directly solving the Schrodinger equation. A fundamental and 
still open question in classical molecular dynamics simulations is how to verify the accuracy computationally, 
i.e., when the solution of the Schrodinger equation is not a computational alternative. 
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The aim of this paper is to derive qualitative error estimates for molecular dynamics and present new 
mathematical methods which could be used also for a more demanding quantitive accuracy estimation, 
without solving the Schrodinger solution. Having molecular dynamics error estimates opens, for instance, 
the possibility of systematically evaluating which density functionals or empirical force fields are good ap- 
proximations and under what conditions the approximation properties hold. Computations with such error 
estimates could also give improved understanding when quantum effects are important and when they are 
not, in particular in cases when the Schrodinger equation is too computational complex to solve. 

The first step to check the accuracy of a molecular dynamics simulation is to know what to compare with. 
Here we compare with the value of any observable g(X), of nuclei positions X, for the time-independent 
Schrodinger eigenvalue equation - H$ = so that the approximation error we study is 

(1.1) / g{X)<P{x,X)*$(x,X)dxdX- lim i / g(X t )dt , 

for a molecular dynamics path X t , with total energy equal to the Schrodinger eigenvalue E. The observable 
can be, for instance, the local potential energy, used in [56] to determine phase-field partial differential equa- 
tions from molecular dynamics simulations, see Figure [T] The time-independent Schrodinger equation has 
a remarkable property of accurately predicting experiments in combination with no unknown data, thereby 
forming the foundation of computational chemistry. However, the drawback is the high dimensional solution 
space for nuclei-electron systems with several particles, restricting numerical solution to small molecules. In 
this paper we study the time-independent setting of the Schrodinger equation as the reference. The pro- 
posed approach has the advantage of avoiding the difficulty of finding the initial data for the time-dependent 
Schrodinger equation. 




Figure 1 . A Lennard- Jones molecular dynamics simulation of a phase transition with pe- 
riodic boundary conditions, from |36| . The left part is solid and the right is liquid. The 
color measures the local potential energy. 



The second step to check the accuracy is to derive error estimates. We have three types of error: time dis- 
cretization error, sampling error and modeling error. The time discretization error comes from approximating 
the differential equation for molecular dynamics with a numerical method, based on replacing time deriva- 
tives with difference quotients for a positive step size At. The sampling error is due to truncating the infinite 
T and using a finite value of T in the integral in ( 1.1 ). The modeling error (also called coarse-graining error) 
originates from eliminating the electrons in the Schrodinger nuclei-electron system and replacing the nuclei 
dynamics with their classical paths; this approximation error was first analyzed by Born and Oppenheimer 
in their seminal paper [2]. 

The time discretization and truncation error components are in some sense simple to handle by comparing 
simulations with different choice of At and T, although it can, of course, be difficult to know that the behavior 
does not change with even smaller At and larger T. The modeling error is more difficult to check since a direct 
approach would require to solve the Schrodinger equation. Currently the Schrodinger partial differential 
equation can only be solved with few particles, therefore it is not an option to solve the Schrodinger equation 
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in general. The reason to use molecular dynamics is precisely in avoiding solution of the Schrodinger equation. 
Consequently the modeling error requires mathematical error analysis. In the literature there seems to be no 
error analysis that is precise, simple and constructive enough so that a molecular dynamics simulation can 
use it to asses the modeling error. Our alternative error analysis presented here is developed with the aim to 
allow the construction of algorithms that estimate the modeling error in molecular dynamics computations. 
Our analysis differs from previous ones by using 

- the time-independent Schrodinger equation as the reference model to compare molecular dynamics 
with, 

- an amplitude function in a WKB-Ansatz that depends only on position coordinates (x,X) (and not 
on momentum coordinates (p, P)) for caustic states, 

- actual solutions of the Schrodinger equation (and not only asymptotic solutions), 

- the theory of Hamilton-Jacobi partial differential equations to derive estimates for the corresponding 
Hamiltonian systems, i.e., the molecular dynamics systems. 

Understanding both the exact Schrodinger model and the molecular dynamics model through Hamiltonian 
systems allows us to obtain bounds for the difference of the solutions by well-established comparison results 
for the solutions of Hamilton-Jacobi equations, by regarding the Schrodinger Hamiltonian and the molecular 
dynamics Hamiltonians as perturbations of each others. The Hamilton-Jacobi theory applied to Hamiltonian 
systems is inspired by the error analysis of symplectic methods for optimal control problems for partial 
differential equations, [30]. The result is that the modeling error can be estimated based on the difference of 
the Hamiltonians, for the molecular dynamics system and the Schrodinger system, along the same solution 
path, see Theorem |5.1| and Section [672] 



2. The Schrodinger and molecular dynamics models 

In deriving the approximation of the solutions to the full Schrodinger equation the heavy particles are 
often treated within classical mechanics, i.e., by defining the evolution of their positions and momenta by 
equations of motions of classical mechanics. Therefore we denote X t : [0, oo) — > R 3N and Pt : [0, oo) — > R 3N 
time-dependent functions of positions and momenta with time derivatives denoted by 

y _ dX t ~ _ d 2 X t 
We denote the Euclidean scalar product on R 3N by 

3N 

X Y = ^2x i Y t . 

i=i 

Furthermore, we use the notation V x^ix, X) = (Vxn/ifi, X), . . . , V^N^fi, X)), and as customary Vj. %p = 

On the other hand, the light particles are treated within the quantum mechanical description and the 
following complex valued bilinear map (•,•): L 2 (R 3n x R 3N ) x L 2 (R 3n x R 3N ) -> L 2 {R 3N ) will be used in 
the subsequent calculations 

(2.1) (0,V> = f ct>(x,X)*iP(x,X)dx. 

The notation ip(x,X) = 0(M~ a ) is also used for complex valued functions, meaning that \ip(x,X)\ = 
0(M~ a ) holds uniformly in x and X. 

The time-independent Schrodinger equation 

(2.2) H(x, X)$(x, X) = E$(x, X) 

models many-body (nuclei-electron) quantum systems and is obtained from minimization of the energy in 
the solution space of wave functions, see [32 [3TJ [TJ [SH [7]. It is an eigenvalue problem for the energy £el 
of the system in the solution space, described by wave functions, $ : R 3n x IR 3 ^ — > C, depending on electron 
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coordinates x — (x 1 , . . . , £™) <E M 3 ", nuclei coordinates X = (X 1 , . . . , X N ) £ R 3N , and the Hamiltonian 
operator H(x, X) 



N 



(2.3) H(x,X)=V(x,X)-hl- 1 Y,A x „. 

n=l 

We assume that a quantum state of the system is fully described by the wave function $ : M. 3n x 
which is an element of the Hilbert space of wave functions with the standard complex valued scalar product 

(($,*))=/" $(x,xyy(x,x)dxdx , 

and the operator T-L is self-adjoint in this Hilbert space. The Hilbert space is then a subset of L 2 (M. 3n x R 3N ) 
with symmetry conditions based on the Pauli exclusion principle for electrons, see [71 122j. 

In computational chemistry the operator V, the electron Hamiltonian, is independent of M and it is 
precisely determined by the sum of the kinetic energy of electrons and the Coulomb interaction between 
nuclei and electrons. We assume that the electron operator V(-,X) is self-adjoint in the subspace with the 



inner product (•, •) of functions in (2.1 ) with fixed X coordinate and acts as a multiplication on functions that 



depend only on X. An essential feature of the partial differential equation (2.2) is the high computational 
complexity of finding the solution in an antisymmetric/symmetric subset of the Sobolev space H 1 (M 3n xM. 3N ). 
The mass of the nuclei, which are much greater than one (electron mass), are the diagonal elements in the 
diagonal matrix M. 

In contrast to the Schrodinger equation, a molecular dynamics model of N nuclei X : [0, T] — > M. 3N , 
with a given potential V p : R 3N — > K, can be computationally studied for large N by solving the ordinary 
differential equations 

(2.4) X t = -V x V p (X t ) , 

in the slow time scale, where the nuclei move 0(1) in unit time. This computational and conceptual 
simplification motivates the study to determine the potential and its implied accuracy compared with the 
the Schrodinger equation, as started already in the 1920's with the Born-Oppenheimer approximation [2]. 
The purpose of our work is to contribute to the current understanding of such derivations by showing 
convergence rates under new assumptions. The precise aim in this paper is to estimate the error 

(2 - 5) j R3N+3 ^(x,xyt>(x,x)dxdx t^Lt ' 9iXt)dt 



for a position dependent observable g(X) of the timc-indepedent Schrodinger equation (2.2 1 approximated 
by the corresponding molecular dynamics observable lim.x->oo T^ 1 J Q T g{X t ) dt, which is computationally 
cheaper to evaluate with several nuclei. The Schrodinger eigenvalue problem may typically have multiple 
eigenvalues and the aim is to find an eigenfunction $ and a molecular dynamics system that can be compared. 
There may be eigenfunctions that we cannot approximate, but with some assumptions on the spectrum of 
V(-,X) the molecular dynamics in fact approximates the observable corresponding to one eigenfunction. 
The main step to relate the Schrodinger wave function and the molecular dynamics solution is the so-called 



zero-order Born-Oppenheimer approximation, where X t solves the classical ab initio molecular dynamics (2.4) 
with the potential V p : M. 3N — >• M determined as an eigenvalue of the electron Hamiltonian V(-, X) for a given 
nuclei position X. That is V P (X) = Xq(X) and 

V(;X)V BO (;X) = \ (X)* BO (;X) , 

for an electron eigenfunction ^>bo{',X) £ L 2 (M. 3n ), for instance, the ground state. The Born-Oppenheimer 
expansion [2] is an approximation of the solution to the time-independent Schrodinger equation which is 
shown in jT51 [TS] to solve the time-independent Schrodinger equation approximately. This expansion, ana- 
lyzed by the methods of multiple scales, pseudo-differential operators and spectral analysis in [TBI H^l [T3"] . 



can be used to study the approximation error (2.5). However, in the literature, e.g., [24] , it is easier to find 
precise statements on the error for the setting of the time-dependent Schrodinger equation, since the stability 
issue is more subtle in the eigenvalue setting. 
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Instead of an asymptotic expansion we use a different method based on a Hamiltonian dynamics formu- 
lation of the time-independent Schrodinger eigenfunction and the stability of the corresponding perturbed 
Hamilton- Jacobi equations viewed as a hitting problem. This approach makes it possible to reduce the error 
propagation on the infinite time interval to finite time excursions from a certain co-dimension one hitting 
set. A motivation for our method is that it forms a sub-step in trying to estimate the approximation error 
using only information available in molecular dynamics simulations. 

The related problem of approximating observables to the time-dependent Schrodinger equation by the 
Born-Oppenheimer expansions is well studied, theoretically in [4J [28] and computationally in |20j using the 
Egorov theorem. The Egorov theorem shows that finite time observables of the time-dependent Schrodinger 
equation are approximated with 0(M~ 1 ) accuracy by the zero-order Born-Oppenheimer dynamics with an 
electron eigenvalue gap. In the special case of a position observable and no electrons (i.e., V = V(X) in 



(2.3)), the Egorov theorem states that 

(2.6) / g[X)*{X,t)**{X,t)dX - f g(X t )*{X Q ,0y$(X ,0)dX 

JR3N J R 3N 

where ^(X, t) is a solution to the time-dependent Schrodinger equation 



< C t M- 



with the Hamiltonian (2.3) and the path X t is the nuclei coordinates for the dynamics with the Hamiltonian 



||X| 2 + y(X). If the initial wave function $(JT, 0) is the eigenfunction in (2.2) the first term in (2.6) reduces 



to the first term in (2.5) and the second term can also become the same in an ergodic limit. However, since 



we do not know that the parameter C t (bounding an integral over (0,i)) is bounded for all time we cannot 



directly conclude an estimate for (2.5) from ( |2.6[ ). 

In our perspective studying the time-independent instead of the time-dependent Schrodinger equation has 
the important differences that 

- the infinite time study of the Born-Oppenheimer dynamics can be reduced to a finite time hitting 
problem, 

- the computational and theoretical problem of specifying initial data for the Schrodinger equation is 
avoided, and 

- computationally cheap evaluation of the position observable g(X) is possible using the time average 
limT->oo Jq g{Xt) dt along the solution path X t . 

In this paper we derive the Born-Oppenheimer approximation from the time-independent Schrodinger 



equation ( 2.2 ) and we establish convergence rates for molecular dynamics approximations to time- independent 
Schrodinger observables under simple assumptions including the so-called caustic points, where the Jacobian 
determinant det J(X t ) = det(dX t /dXo) of the Eulerian-Lagrangian transformation of X-paths vanish. As 
mentioned previously, the main new analytical idea is an interpretation of the time-independent Schrodinger 



equation (2.2) as a Hamiltonian system and the subsequent analysis of the approximations by comparing 
Hamiltonians. This analysis employs the theory of Hamilton- Jacobi partial differential equations. The prob- 
lematic infinite-time evolution of perturbations in the dynamics is solved by viewing it as a finite-time hitting 
problem for the Hamilton-Jacobi equation, with a particular hitting set. In contrast to the traditional rigor- 
ous and formal asymptotic expansions we analyze the transport equation as a time-dependent Schrodinger 
equation. 

The main inspiration for this paper are works [37J G3 IS] and the semi-classical WKB analysis in [25] : 
the works (27J [6l [5] derive the time-dependent Schrodinger dynamics of an x-system, ity = Hi^>, from 
the time- independent Schrodinger equation (with the Hamiltonian Tli(x) + el-L{x, X)) by a classical limit 
for the environment variable A, as the coupling parameter e vanishes and the mass M tends to infinity; 
in particular [HI [S] show that the time derivative enters through the coupling of ^> with the classical 
velocity. Here we refine the use of characteristics to study classical ab initio molecular dynamics where 
the coupling does not vanish, and we establish error estimates for Born-Oppenheimer approximations of 
Schrodinger observables. The small scale, introduced by the perturbation 



of the potential V, is identified in a modified WKB eikonal equation and analyzed through the corresponding 
transport equation as a time-dependent Schrodingcr equation along the eikonal characteristics. This modified 
WKB formulation reduces to the standard semi-classical approximation, see [25 , in the case of the potential 
function V = V(X) G K, depending only on nuclei coordinates, but becomes different in the case of operator- 
valued potentials studied here. The global analysis of WKB functions was initiated by Maslov in the 
I960', [2"5] , and lead to the subject Geometry of Quantization, relating global classical paths to eigenfunctions 
of the Schrodinger equation, see 10J . The analysis presented in this paper is based on a Hamiltonian system 
interpretation of the time-independent Schrodinger equation. Stability of the corresponding Hamilton- Jacobi 
equation, bypasses the usual separation of nuclei and electron wave functions in the time-dependent self- 
consistent field equations, [51 [2"31 155] . 

Theorem |5.f | demonstrates that observables from the zero-order Born-Oppenheimer dynamics approxi- 
mate observables for the Schrodinger eigenvalue problem with the error of order 0(M~ 1+S ), for any S > 0, 
assuming that the electron eigenvalues satisfy a spectral gap condition. The result is based on the Hamil- 



tonian (2.3) with any potential V that is smooth in X, e.g., a regularized version of the Coulomb potential. 
The derivation does not assume that the nuclei are supported on small domains; in contrast derivations 
based on the time-dependent self-consistent field equations require nuclei to be supported on small domains. 
The reason that the small support is not needed here comes from the combination of the characteristics 
and sampling from an equilibrium density. In other words, the nuclei paths behave classically although 
they may not be supported on small domains. Section [4] shows that caustics couple the WKB modes, as is 
well-known from geometric optics, see [18; 25!, and generate non-orthogonal WKB modes that are coupled 
in the Schrodinger density. On the other hand, with a spectral gap and without caustics the Schrodinger 
density is asymptotically decoupled into a simple sum of individual WKB densities. Section [7] constructs a 



WKB-Fourier integral Schrodinger solution for caustic states. Section 5.2 relates the approximation results 
to the accuracy of symplectic numerical methods for molecular dynamics. 

A unique property of the time-independent Schrodingcr equation we use is the interpretation that the 
dynamics X t g M. 3N can return to a co-dimension one surface I which then can reduce the dynamics to 
a hitting time problem with finite-time excursions from /. Another advantage of viewing the molecular 
dynamics as an approximation of the eigenvalue problem is that stochastic perturbations of the electron 
ground state can be interpreted as a Gibbs distribution of degenerate nuclei-electron eigenstates of the 



Schrodinger eigenvalue problem (2.2), see 33 . The time-independent eigenvalue setting also avoids the issue 
on "wave function collapse" to an eigenstate, present in the time-dependent Schrodinger equation. 

We believe that these ideas can be further developed to better understanding of molecular dynamics 
simulations. For example, it would be desirable to have more precise conditions on the data (i.e. molec- 
ular dynamics initial data and potential V) instead of our implicit assumption on finite hitting time and 



convergence of the Born-Oppenheimer power series approximation in Lemma |6.2| 

3. A TIME-INDEPENDENT SCHRODINGER WKB-SOLUTION 

3.1. Exact Schrodinger dynamics. For the sake of simplicity we assume that all nuclei have the same 
mass. If this is not the case, we can introduce new coordinates Ml X k = Mu' 2 X k , which transform the 
Hamiltonian to the form we want V{x,M^ 2 M~ l t 2 X) — (2M%)~ 1 J2k=i ^-x k ■ The singular perturbation 
— (2M) _1 ^2 k A X k of the potential V introduces an additional small scale M~ x l 2 of high frequency oscilla- 



tions, as shown by a WKB-expansion, see [2S1 El HE1 26 . We shall construct solutions to (2.2) in such a 
WKB-form 

(3.i) $(.t,a) = <f>( x ,xy Ml/2e ( x \ 

where the amplitude function <f> : K 3n x M. 3N — > C is complex valued, the phase 9 : R 3N — > K is real 
valued, and the factor M 1 ' 2 is introduced in order to have well-defined limits of <ft and 9 as M — » oo. 



Note that it is trivially always possible to find funtions <f> and 9 satisfying (3.1 ), even in the sense of a true 
equality. Of course, the ansatz only makes sense if (f> and 9 do not have strong oscillations for large M. 
The standard WKB-construction, j^SJITB], is based on a series expansion in powers of M 1 / 2 which solves 
the Schrodinger equation with arbitrary high accuracy. Instead of an asymptotic solution, we introduce an 
actual solution based on a time-dependent Schrodinger transport equation. This transport equation reduces 
to the formulation in 25 for the case of a potential function V = V(X) £ R, depending only on nuclei 
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coordinates X £ K 3Ar , and modifies it for the case of a self-adjoint potential operator V(-, X) on the electron 
space L 2 (M 3 ™) which is the primary focus of our work here. In Sections [4] and [7] we use a linear combination 
of WKB-eigensolutions, but first we study the simplest case of a single WKB-eigensolution as motivated by 
the following subsection. 



3.1.1. Molecular dynamics from a piecewise constant electron operator on a simplex mesh. The purpose of 
this section is to convey a first formal understanding of the relation between ab initio molecular dynamics 
X t = — Vx^o(Xt) and the Schrodinger eigenvalue problem (2.2| and motivate the WKB ansatz ( 3.1 1. In 
subsequent sections we will describe precise analysis of error estimates for the WKB-method. The idea 
behind this first study is to approximate the electron operator V by a finite dimensional matrix V h , which 
is piecewise constant on a simplex mesh in the variable X, with the mesh size h. Furthermore, we introduce 
the change of variables 

./ 

3=0 

based on the piecewise constant electron eigenvalues and eigenvectors V hx Sj = ^'j^j, (^j^j) = 1> j = 
0, ... J, normalized and ordered with respect to increasing eigenvalues. Then the Schrodinger equation (2.2) 
becomes 

with the notation = y\ , so that on each simplex 

~2~M Axtft + = E<Pj ' 
which by separation of variables, for each j = 0, 1, 2, . . . , J, implies 

AM 1/2 P j -X 

uyjr- )c 

Pi 

for any P- 3 G C 3N that satisfies the eikonal equation 



(3.2) ^ = ^o(i«y 



\p3 ■ P j + X 1 * = E , 

for any a(P- 7 ) £ C, if all components of pi are non zero. If P J k — we have a(P : ') = Yi{k ■ p j =o} i^-k^k + Bf.) 
for any Ak G C, B^ 6 C, since e ±-lMl,2pl k x k = \ m this case . The solution to (2.2 1, and its normal 
derivative are continuous at the interfaces of the simplices. On the intersection of the faces the normal 
derivative is not defined but this set is of measure zero and thus negligible as seen from the iP(]R 3Ar ) 
solution concept of (2.2). 

We investigate a simpler, one-dimensional case, l€K, first. Then the solution ip simplifies to 

iM 1/2 P j -X . 1 -iM 1/2 P j -X 

for dj, bj, Pi e C and (P-?) 2 /2 + Xj — E . The continuity conditions 

lim = lim 

[ ' ' lim d x <5>{X) = lim d x §{ x ) 

X— fXn+ X—>Xq — 

hold for any Xq G R, in particular, at the interval boundary where for Xq = 

(3.4) 3 

. lim ^ a x $(X) = iM 1 / 2 V(a J± P^ ± - b J± P^ j± ) . 



E' 

j 
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It is clear that given a_ and 6_ we can determine a+ and &+ so that (3.3 1 holds. In order to prepare for the 



multi-dimensional case it is convenient to consider each incoming wave a_ and b + separately: the incoming 
a_ wave is split into a refracted a + and reflected 6_ wave 

(3.5) "■■ ( " - / " = £K-+^-+4 + bj-Vj-PL) 

j j 

and similarly the incoming b + wave is split into a refracted 6_ wave and a reflected a + wave, see Figure [2] 
The jump conditions at the different interfaces are coupled by the oscillatory functions e ±tM 1 pl x . The 
global construction of tp and ^ in one dimension follows by marching in the positive A-direction to successive 
intervals, creating in each interval both a e rM 1 pJ ' x ^j and a e~ zM ' pJ ' x tyj wave. 

In general each interface condition (3.4 1 also couples all eigenvectors However, we shall see that if 
M is large, V smooth and there is a spectral gap Ai — Ao > c > then, in the limit of the simplex size h 
tending to zero, there is an asymptotically uncoupled WKB-solution $(x,X) — <p(x, X)e lM 1 ( x \ where 
9 : R 3 N — » R, 4> : ]R 3n x M. 3N —> C. Under these assumptions the Born-Oppenheimer approximation in 
Lemma 6.2 shows that <f> is asymptotically parallel, in L 2 (dx), to the electron eigenfunction as M —> oo. 
The gradient V X 6{X) = P° is obtained from the differential 0(A) = 0(A o )+V ' x 9(X Q )-(X-X )+o{\X-X Q \). 

In the case of electron eigenvalue crossing, i.e., Ai(A) = Xq(X) for some X, or so called avoided crossings 
(meaning that the eigenvalue gap c -C 1 is small and dependent on M) , a refraction will, in general, include 

all components aje lM p3 x ^j, j = 1, . . . , J and consequently the Born-Oppenheimer approximation fails. 

The construction of a solution to the Schrodinger equation with a piecewise constant potential is more 
involved in the multi-dimensional case for two reasons: each reflection at an interface generates, in general, 
an additional path in a new direction, so that many paths are needed. Furthermore, the construction of a 



solution to the eikonal equation is more complicated since the jump condition ( 3.4 ) implies that the tangential 
component P/ of P 3 must be continuous across a simplex face and only the normal component P-{ = P J 



PI 



may have a jump. In multi-dimensional cases it is still possible to construct a solution of the form (3.2) by 
following the characteristic paths X t = P 3 (X t ) and using the jump conditions (3.4): when the path X t hits 



a simplex face, the tangential part Pl of P 3 is continuous and the normal component P^ of P 3 may jump. 
At a simplex face the new value of the P^ is determined by (P 3 t ■ P 3 t + Pj! ■ -P/)/2 + A^ = E. Analogously 
to the one dimensional case we treat the pair e tM± 2 ( p t+ p n)- x and e lMl/2 ( p * J ~ P ")' x together. However, 
each collision with e tM 1 ( p t+ p ™)' x on an interface now creates a reflected wave in another direction, in 
particular, e lM ' ( p t ,_p n)' x v[/ J ; , and we get many paths to follow. Therefore each mode e lM ' pl ' x follows 
its characteristic X t , where X t — P 3 , through the simplex to the adjacent simplicial faces, which the 
characteristic pass through when they leave the simplex, and at these outflow faces a reflected mode is 
created and a refracted mode continues into the adjacent simplices, see Figure [2j In this way we can 



formally construct a solution of the form ^2 P j a(P 3 )e^ ll/2pl ' Xl >f?j to the Schrodinger equation (2.2), with 
possibly several different characteristic paths in each simplex. 




Figure 2. The value of P 3 is constructed by following the characteristic paths X t (the blue 
and green curves), based on X t — P 3 , with a reflection-refraction at each simplex face (left) 
following the path through simplices (middle) and each simplex may have several P 3 (right). 



In conclusion, the piecewise constant electron operator shows that the solution to the Schrodinger equation 
(2.2) is composed of a linear combination of highly oscillatory function modes a 3 'e* 1 3p '' X $j based on the 
electron eigenvectors \Pj and eigenvalues Xj, where P 3 satisfies the eikonal equation P 3 ■ P 3 /2 + \j(X) = E. 



These modes can be followed be characteristics X = P^ from simplex to simplex. In this paper we show that 
observables based on the related WKB Schrodinger solutions can be approximated by molecular dynamics 
time averages, when there is a spectral gap around Aq. 



3.1.2. A first WKB-solution. The WKB-solution satisfies the Schrodinger equation (2.2) provided that 

= (H - E)</> e iMl/2(> W 



(3.6) 



V-E)(j> 



1 

2M' 



JM 1/2 9(X) 



We shall see that only eigensolutions $ that correspond to dynamics without caustics correspond to such 
a single WKB-mode, as for instance when the eigenvalue E is inside an electron eigenvalue gap. Solutions 
in the presence of caustics use a Fourier integral of such WKB-modes, and we treat this case in detail in 
Section [7] To understand the behavior of 6 >, we multiply (3.6) by <f>* e ~ tMl/2 Q( x ) anc [ integrate over M 3n . 
Similarly we take the complex conjugate of (3.6), and multiply by (f>e lM 1 9 ^ x ^ and integrate over M 3 ". By 
adding these two expressions we obtain 



= 2{-\V x t 



-E) 



(3.7) 



MV2 v 



(Vxf V x 0,<t>)) 



2M 



2MV2 



({ct>, 4>) - {&<!>)) A 



=2iIm(0,Vx0-V x e) =0 

The purpose of the phase function 9 is to generate an accurate approximation in the limit as M — > 00. A 
possible and natural definition of 9 would be the formal limit of (3.7) as M — > 00, which is the Hamilton- 
Jacobi equation, also called the eikonal equation 

(3.8) l\V x 9\ 2 =E-V , 



where the function Vq 
(3.9) 



V 



{<j>,V<j>) 



The solution to the Hamilton-Jacobi eikonal equation can be constructed from the associated Hamiltonian 
system 

X t = P t 

(3.10) 

Pt = -VxV (X t ) 

through the characteristics path (X t ,P t ) satisfying Vj9(Xj) =: Pt- The amplitude function can be 
determined by requiring the ansatz (3.6) to be a solution, which gives 

(H-E)0e tMl/2() W 







'At 



V -E), 



1 

~2M 

so that by using (3.8) we have 

1 



A x cf> + (V-Vo)(i>- 



M 1 / 2 



(V X (t>-V x + \<f>A 



At 



a iM 1/2 8(X) 



A x ^+(V-V )cj>- 



:(V 



A< 



V; 



-(f) A 



At 



= 0. 



The usual method for determining <fi from this so-called transport equation uses an asymptotic expansion 
~ J2k=o M~ k/2 <f>k, see [HUH] and the be ginning of Section^ An alternative is to write it as a Schrodinger 
equation, similar to work in |25| : we apply the characteristics in ( |3. 10 1 to write 

j t <KXt) = Va</> ■ x t = v x <t> ■ v x e , 



and define the weight function G by 
(3.11) 



2 &xO(X t ) , 



and the variable ip t '■= (/>(X t )Gt- We use the notation <f>{X) instead of the more precise </>(•, X), so that e.g. 
ifj t — ipt{ x ) = <t>( x >X t )G t . Then the transport equation becomes a Schrddinger equation 



(3.12) 



iM- 1 ^ = (V - V )ip t 



Gt_ 
2M" 



A 



x 



In conclusion, equations (3.8)-(3.12) determine the WKB-ansatz (3.1| to be a solution to the Schrddinger 
equation ( 2.2 1. 

Theorem 3.1. Assume the Hamilton- Jacobi equation, with the corresponding Hamiltonian, 



H S (X,P) :=i|P| 2 



mX),V(X)iP(X)) 

mx)^j(x)) 

S v ' 

=:V (X) 



-E = 0, 



based on the primal variable X and the dual variable P = P(X) = VxJ(X), has a smooth solution 6{X), 
then 9 generates a solution to the time-independent Schrddinger equation (jH — E)<& = 0, in the sense that 



$(X t ,x) = G- 1 (X t )*l>(x,X t )e 



iM 1/2 6(X t ) 



solves the equation \2.2\, where ip{Xt) := ipt satisfies the transport equation (3.12 I and 



G(X t 
d 



G, 



1 



dt io gGt = -A x e { x t ) , 



(X t ,Pt) solves the Hamiltonian system (3.10) corresponding to Hg- 

It is well know that Hamilton- Jacobi equations in general do not have smooth solutions, due to X-paths 
that collide, as seen by (7.25) generating blow up in dxxO(X). However if the domain is small enough, the 
data on the boundary is smooth and Vo is smooth, then the characteristics generate a smooth solution, see 
Ref. [12]. In Section 7.2.5 we describe Maslov's method to find a global solution by patching together local 
solutions. 

Note that the nuclei density, using G, can be written 
(3.13) p 



Jr 



• dX 



l R 3N \V) vi »>■<*■ J R3N (tp, i>)G 2 dX 

and since each time t determines a unique point (X t ,P t ) = (X t ,Vx9(X t )) in the phase space the functions 
G and ip are well defined. 

3.1.3. Liouville's Formula. In this section we verify Liouville's formula 

d(X ) 



(3.14) 



G\ 



det 



d(X t ) 



given in [25]. The characteristic X t = P(X t ) implies f t J(X t ) = V X PJ(X), where J{X) tj = dX l t /dX J 
denotes the first variation with respect to perturbations of the initial data. The logarithmic derivative then 
satisfies d/dt( log J(X)) = dxjP l {X t ) = dxixoO{X) which implies that log J[Xt) is symmetric and shows 
that ( |3.14[ ) holds 

divP = Tr V X P = ^-Tr log J(X) = ^ log det J(X) . 
dt dt 

The last step uses that J(X) can be diagonalized by an orthogonal transformation and that the trace is 
invariant under orthogonal transformations. 

in 



3.1.4. Data for the Hamiltonian system. For the energy E chosen larger than the potential energy, that is 
such that E > Vq, the Hamiltonian system (3.10) yields a solution (X,P) : [0,T] — > U x R 3N to the cikonal 
equation (3.8) locally in a neighborhood U C R :iN , for regular compatible data (Xq,Po) given on a 3iV — 1 



dimensional " inflow" -domain I C.U. Typically, the domain / and the data (Xq,Pq)\i are not given (except 
that its total energy is E), unless it is really an inflow domain and characteristic paths do not return to / as 
in a scattering problem. If paths leaving from / return to /, there is an additional compatibility of data on 
I: assume X 6 / and X t £ I, then the values P t are determined from P ; continuing the path to subsequent 
hitting points X tj <G I, j = 1, 2, . . . determines Pt j from Po. The characteristic path (Xt, Pt), t > 0, generates 
a manifold in the phase space (X, P), which is smooth under our assumptions. This manifold is in general 
only locally of the form (X, P(X)), but in the case of no caustics it is globally of this form and then there 
is a phase function X H> 9(X) such that P(X) = V x@(X) globally. In Section [7] we study phase space 
manifolds with caustics. 

Remark 3.2. The integrating factor G and its derivative dx^G can be determined from (P, dx*P, 9x^x3 P) 



along the characteristics by the following characteristic equations obtained from (3.8 1 by differentiation with 
respect to X 



d 
dt 



d X rP k 



(3.15) 



-j,dx^xiP k 
at 



J2P 3 d X lXrP k =J2 P ' d XrX«P 3 

3 3 

J2dxrP j d x *P j -dxrxxVo, 
j 

Yl p3 dxixrx«P k + pjd xrx*x«P j 

3 3 

Y d xrP 3 dx* Xq P 3 - J2 d xrx«P j d x *P 3 - d XrX *x«V , 



and similarly dx^xiG can be determined from (P, dxiP, dx^xiP, dxixix k P)- 

3.2. Born-Oppenheimer dynamics. The Born-Oppenheimer approximation leads to the standard for- 
mulation of ab initio molecular dynamics, in the micro-canonical ensemble with the constant number of 
particles, volume and energy, for the nuclei positions X — X-qq, 



(3.16) 



X t =Pt, 

Pt = -VxAopft) , 



by using that the electrons are in the eigenstate rp — ^bo with eigenvalue Ao to V, in L 2 (dx) for fixed X, 
i.e., V(X)t>BO = A (X)*bo- The corresponding Hamiltonian is H BO (X,P) := |P| 2 /2 + \q(X) with the 
eikonal equation 



(3.17) 

3.3. Equations for the density. We note that 

4> = G~ 1 4> 



-\V x 0bo(X)\ 2 + X (X) = E. 



1/2 



shows that G and tp determine the density 
(3.18) p s - P 



J(i>,iP)\G\- 2 dX 



defined in (3.13). Using the Born-Oppenheimer approximation in Lemma 6.2 we have (ip, ip) — 1 + 0(AI^ 1 ) 
in the case of a sp ectral gap. Therefore the weight function |G| -2 approximates the density and we know 
that |G| -2 is determined by the phase function 9. 

li 



from Theorem 



3.1 



The Born-Oppenheimer dynamics generates an approximate solution ^boG B qC iM ' 9bo which yields the 
density 

(3.19) p BO = \Gbo\~ 2 , 
where 

j t \og\G BO \- 2 = -A x 6 BO (X) . 
This representation can also be obtained from the conservation of mass 

(3.20) = div(pBoVx^Bo) 
implying 

(3.21) ^rPBo(^t) = VxPBo(Xt) ■ X t = -pBo(X t ) divV x 6 BO , 

at 

with the solution 

(3.22) PBo(Xt) ( 



\G BO (X t )\ 2 

where C is a positive constant for each characteristic. Note that the derivation of this classical density does 
not need a corresponding WKB equation but uses only the conservation of mass that holds for classical paths 
satisfying a Hamiltonian system. The classical density corresponds precisely to the Eulerian-Lagrangian 
change of coordinates | G t \ 2 / \ G \ 2 = det (dX t /dX ) in (jglg) . 



3.4. Construction of the solution operator. The WKB Ansatz (3.1) is meaningful when tp does not 



include the full small scale. In Lemma |6.2| we present conditions for ijj to be smooth. 

To construct the solution operator it is convenient to include a non interacting particle in the system, i.e., 
a particle without charge, and assume that this particle moves with a constant, high speed dX\/dt = ^> 1 
(or equivalently with the unit speed and a large mass). Such a non interacting particle does not affect the 
other particles. The additional new coordinate X\ is helpful in order to simply relate the time-coordinate t 
and X\. We add the corresponding kinetic energy (P/) 2 /2 to E in order not to change the original problem 
(2.2) and write the equation (3.12) in the fast time scale r = M x l 2 t 

= (V - W " 27tf G E Axi (G ^ } • 
j 

Furthermore, we change to the coordinates 

(t,X*) := (t,X%,X$,X 2 ,...,X n ) G [0,oo) x I, instead of {X\ X 2 , . . . , X N ) e R m , 
where X 1 = (X{,X J 2 ,X J 3 ) e M 3 . Hence we obtain 

(3-23) ii, + ^rp^ - (V - V )^ - ^ A xi (G~ V) =: Vty , 



using the notation w = dw/dr in this section. In Section 6.1 we show that the left hand side can be reduced 
to itp as Pi — > oo, by choosing special initial data. Note also that G is independent of the first component 
in X 1 . We see that the operator 



V := G-'VG = G-\V 7 )G-JL£ A xi 



2M 

=v-v 3 

is symmetric on J L 2 (M 3n+3Ar - 1 ). Assume now the data (A ,P ,^o) for X e E 3Ar_1 is (LZ^^-periodic, 
then also (X T ,P T ,Z T ) is ( J LZ) 3Ar " 1 -periodic, for Z t = 6(X t ) and P t = V x 0{X t ). To simplify the notation 
for such periodic functions, define the periodic circle 

T :=R/(LZ). 

We seek a solution $ of (2.2) which is (LZ) 3 ( n+Ar ^ 1 -periodic in the (x, X*)-variable. The Schrodinger 
operator V(-, X T ) has, for each r, real eigenvalues {A m (r)} with a complete set of eigenvectors {( m (x, X*, r)} 
orthogonal in the space of x-anti-symmetric functions in L 2 (T 3n+3JV_1 ), see pQ. The proof uses that the 
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operator V T + 7/ generates a compact solution operator in the Hilbert space of x-anti-symmetric func- 
tions in L 2 (T 3n+3Ar_1 ), for the constant 7 £ (0, 00) chosen sufficiently large. The discrete spectrum and 
the compactness comes from Fredholm theory for compact operators and the fact that the bilinear form 
Jt3("+«)- 1 v VtW + 'yvw dx dX* is continuous and coercive on H 1 (T 3 ^ n+N ^~ 1 ), see [T^]. We see that V has the 
same eigenvalues {A m (r)} and the eigenvectors {G T C m (r)}, orthogonal in the weighted L 2 -scalar product 



(v,w) G-'-dX* . 

The construction and analysis of the solution operator continues in Section [6 . 1 1 based on the spectrum 



Remark 3.3 (Boundary conditions). The eigenvalue problem (2.2) makes sense not only in the periodic 
setting but also with alternative boundary conditions from interaction with an external environment, e.g., 
for scattering problems. 



4. Computation of observables 
Suppose the goal is to compute a real-valued observable 



[ (<$>,A$)dX 

Jf3N 



for a given bounded linear multiplication operator A = A(X) on L 2 (T 3N ) and a solution $ = J2k 4 l k etMl/29k 



of J2J2J). We have 
(4.1) 



= Em Jt3« Ae^^W^W)^, dX . 



The integrand is oscillatory for k =^ I, hence critical points (or near critical points) of the phase difference 
give the main contribution. The stationary phase method, see [TOl US] and Section [9j shows that these 
integrals are small, bounded by 0(A/~ 3Ar / 4 ), in the case when the phase difference has non degenerate 
critical points, or no critical point, and the functions A(fa,(j>i) and 9i are sufficiently smooth. A critical 
point X c £ K 3JV satisfies V x&i{X c ) — V x9k{X c ) = 0, which means that the two different paths, generated 
by 9i and 9k, passing through X — X c also have the same momentum P at this point. That the critical 
point is degenerate means that the Hessian matrix dx^xj ($fc — Qi)(X c ) is singular (or asymptotically singular 
for M — ► 00 as for avoided crossings when the electron eigenvalues have a vanishing spectral gap depending 
on M). Therefore caustics, crossing or avoided crossing electron eigenvalues may generate coupling between 
the WKB terms. On the other hand, without such coupling the density of a linear combination of WKB 
terms separates asymptotically to a sum of densities of the individual WKB terms 



(4.2) 



k 

[ ($, A$)dX = V / A (fa, fa) dX + OiM- 1 ) 
Jt 3N Jt 3N ~* — v — ' 



in the case of multiple eigenstates, k > 1, and 

/ (<P,A^)dX= [ A(fa,fa)dX 
Jt 3N Jt 3N 

for a single eigenstate. In the next section we will study molecular dynamics approximations of a single state 

(4.3) / A(fa,fa)dX= [ A(X)p k {X)dX. 

7t 3N 7t 3N 

In the presence of a caustic, the WKB terms can be asymptotically non orthogonal, since their coefficients 
and phases typically are not smooth enough to allow the integration by parts to gain powers of M^ 1 / 2 . 
Non-orthogonal WKB functions tell how the caustic couples the WKB modes. 

Regarding the inflow density Pk\r there are two situations: either the characteristics return often to the 
inflow domain or not. If they do not return we have a scattering problem and it is reasonable to define the 
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inflow-density pt\j as an initial condition. If characteristics return, the dynamics can be used to estimate 
the return-density pk\ t as follows: Assume that the following limits exist 



(4.4) lim - / A{X t )dt= / A(X) Pk (X)dX 



T 



J3N 



which bypasses the need to find pt\ I and the quadrature in the number of characteristics. A way to think 
about this limit is to sample the return points X t £ I and from these samples construct an empirical return- 
density, converging to pk \ I as the number of return iterations tends to infinity. We shall use this perspective 



to view the cikonal equation (3.8) as a hitting problem on /, with hitting times t (i.e., return times). The 
property having p\j constant as a function of A is called ergodicity, which we will use. We could allow 
the density p\j to depend on the initial position Xq and momentum Pq, but then our observables need to 
conditional expected values. An example of a hitting surface is the co-dimension one surface where the first 
component X\\ in X\ — (An, A12, A13) is equal to its initial value An(0). The dynamics does not always 
have such a hitting surface: for instance if all particles are close initially and then are scattered away from 
each other, as in an explosion, no co-dimension one hitting surface exists. 



5. Molecular dynamics approximation of Schrodinger observables 



c , Acj)k) dX has the main ingredients: 



A numerical computation of an approximation to ^ fc J T3N 

(1) to approximate the exact characteristics by molecular dynamics characteristics (3.10), 

(2) to discretize the molecular dynamics equations, and 

(3a) if p\j is an inflow-density, to introduce quadrature in the number of characteristics, or 



(3b) if p\ is a return-density, to replace the ensemble average by a time average using the property (4.4) 

This section presents a derivation of the approximation error in the step (1) in the case of a return density and 
comments on the time-discretization of step (2) treated in Section |5. 2 1 The third and fourth discretization 
steps, which are not described here, are studied, for instance, in 



5.1. The Born-Oppenheimer approximation error. This section states our main result of molecu- 
lar dynamics approximating Schrodinger observables. We formulate it using the assumption of the Born- 
Oppenheimer property 



(5.1) 



IIV* - *Bo(*t)IU»(«te) = 0(M 



-1/2N 



uniformly in t. 



This assumption is then proved in Lemma [6. 2| based on a setting with a spectral gap. 

The spectral gap condition. The electron eigenvalues {A*,} satisfy, for some positive c, the spectral gap 
condition 



(5.2) 



inf |A fe (y)-A (Y)|>c, 



where D :— {Ag(i) 1 1 > 0} U {Abo(^) 1 1 > 0} is the set of all nuclei positions obtained from the Schrodinger 
characteristics X — As in Theorem 3.1 and from the Born-Oppenheimer dynamics X = Abo in (3.16), for 
all considered initial data. 

Theorem 5.1. Assume that the phase functions 9s and #bo ar £ smooth solutions to the eikonal equations 



(3.8) and (3.17) and that the Born-Oppenheimer property (5.1) holds, then the zero-order Born-Oppenheimer 



dynamics (3.16), assumed to have the ergodic limit (4.4) and bounded hitting times r in (6.11), (6.14) and 



(6.18), ap proxi mates time-independent Schrodinger observables, generated by Theorem 3.1 or the caustic case 
in Section 7.2 with error bounded by 0(M~ 1+S ) 



(5.3) 



J3N 



g(X) PBO (X)dX 



J3N 



g(X)p s (X) dX + 0{M- 1+S ) , for any 5 > 0. 



The proof is given in Sections [6] and |7.2| 
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5.2. Why do symplectic numerical simulations of molecular dynamics work? The derivation of the 
approximation error for the Born-Oppenheimer dynamics, in Theorem |5.1[ also allows to study perturbed 
systems. For instance, the perturbed Born-Oppenheimer dynamics 

X t = P t +VpH e (X t ,Pt) 

Pt = -V x \ {X t ) - V x H e (X t ,P t ) , 
generated from a perturbed Hamiltonian Hbo(X, P) + H e (X,P) = E, with the perturbation satisfying 
(5.4) ||-H" e ||i- < e for some e G (0, oo) 



yields through (6.13) and (6.19) an additional error term 0(e) to the approximation of observables in (5.3). 
So called symplectic numerical methods are precisely those that can be written as perturbed Hamiltonian 
systems, see |30j . and consequently we have a method to precisely analyze their numerical error by combin 



ing an explicit construction of H e with the stability condition (5.4) to obtain 0((M 1 + e) 1 5 ) accurate 



approximations, provided the corresponding phase function has bounded second difference quotients. The 
popular Stormer-Verlet method is symplectic and the positions X coincides with those of the symplectic 
Euler method, for which H e is explicitly constructed in [30 with e proportional to the time step. The 
construction in |30j is not using the modified equation and formal asymptotics, instead a piecewise linear 
extension of the solution generates H e . 

6. Analysis of the molecular dynamics approximation 

Before we proceed with the analysis of the approximation error we motivate our results by a significantly 
simpler case of a system without electrons. We use the densities (3.18) and ( 3.19| ) and we show heuristically 



how the characteristics can be used to estimate the difference ps — Pbo, leading to 0(M x ) accurate Born- 
Oppenheimer approximations of Schrodinger observables 

g(X) Ps (X) dX = J g(X)p BO (X)dX + 0(M- 1 ). 

<*,*> 

In the special case of no electrons, the dynamics of X does not depend on ip and therefore X^o = X§ = X 
and consequently Gbo = @S ■ The difference ip$ — "0BO can De understood from iterative approximations of 
d3~12 ) 

(6-1) j^72^+i ~ { y~ = ^GAx(G-Vfc) 

with ipQ = 0. Then V'bo = V'l is the Born-Oppenheimer approximation and formally we have the iterations 
approaching the full Schrodinger solution ipf. ips as k — > oo. 

In the special case of no electrons, there holds V = Vo, thus the transport equation iipi = has constant 
solutions. We let tpi = 1 and then i/> 2 — ipi is imaginary with its absolute value bounded by 0(M~ 1 / 2 ). We 



write the iterations of ipk by integrating (6.1 ) as the linear mapping 

k 

^ fe+1 = l+iM-^S^u) = Y^i l M' l/2 S\^i) , 

1=0 

which formally shows that 

|Vs| 2 = IV'il 2 + 2Re (tfa - Vi,^i> + |Vs - ^Ail 2 = 1 + 0(M~ l ) . 
Consequently this special Born-Oppenheimer density satisfies 

(6.2) Pbo = Gg 2 (^s>s) +0(M- 1 ) , 

since Gbo = Gs and X do not depend on ip. 

In the general case with electrons and a spectral gap, we show in Lemma 6.2 that there is a solution 
satisfying 

(6.3) Us-^Bo\\m dX ) ^o(Ar^ 2 ), 
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for the electron cigenfunction $bo, satisfying 

V(;X)V BO (;X) = A (X)* BO (;X) 

and the eigenvalue \q(X) G M with a (fixed) nuclei position X. Then the state ipi equal to a constant, in 
the case of no electrons, corresponds to the electron eigenfunction ^bo in the case with electrons present. 
In the general case the X dynamics for the Schrodinger and the Born-Oppenheimer dynamics are not the 
same, but we will show that (6.3) implies that the Hamiltonians H$ and Pbo are 0(M~ 1 ) close. Using 
stability of Hamilton-Jacobi equations, the phase functions 9$ and #bo arc then also close in the maximum 
norm, which, combined with an assumption of smooth phase functions, show that |Gs — Geo I = 0(M~ 1+S ) 

also shows that l(V'S)V's) — 1| = (^(AI^ 1 ) and consequently the density bound 



6.2 



for any S > 0. Lemma 

|ps — Pbo| = 0(M~ 1+S ) holds. To obtain the estimate (6.3) the important new property, compared to no 
electrons, is to use oscillatory cancellation in directions orthogonal to ^bo- 

6.1. Continuation of the construction of t he s olution operator. This section continues the construc- 



tion of the solution operator started in Section 3.4 Assume for a moment that V is independent of r. Then 



the solution to (3.23) can be written as a linear combination of the two exponentials 

ae iTA + + be lTA - 



where the two characteristic roots are the operators 

A ± = (Pi) 2 (-1 ± (1 - 2(P 1 1 )- 2 V) 
We see that e lTA - is a highly oscillatory solution on the fast r-scale with 



1/2 



lim t— r 



2Id, 



while 
(6.4) 

Therefore we chose initial data 
(6.5) 



lim A+ 



-V. 



oo determines the solution by the 



#| T =o = -A + ip\ T =a 

to have b — 0, which eliminates the fast scale, and the limit 
Schrodinger equation 

iip = Vtp . 

The next section presents an analogous construction for the slowly, in r, varying operator V. 



6.1.1. Spectral decomposition. Write (3.23) as the first order system 

iii = -2(P 1 1 ) 2 (VV -tt) , 



which for ip :— (ip,Tr) takes the form 



-1 
2(P 1 1 ) 2 V -2(pi) 2 

where the eigenvalues A± , right eigenvectors Q± and left eigenvectors Q^ 1 of the real "matrix" operator B 
are 

1/2N 



A± := (PIY -Id ± Id - 2(P/)- 2 V 



Id 

-A+ 



Q; 1 := (Id-A+Al 1 



Q- 



-Al 
Id 



Id 

a: 1 



q: 1 
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(id-A + (A_r 1 )- 1 



A+ 
Id 



We see that limpi^^ A+ = —V and lim Pi i^ co (P 1 1 ) _2 A_ = — 2Id. The important property here is that 
the left eigenvector limit limpi^^ Q^ 1 = (Id, 0) is constant, independent of r, which implies that the Q+ 
component Q^ 1 ^ = tp decouples. We obtain in the limit — > oo the time-dependent Schrodinger equation 

or dr 

= -A+(r)Q; 1 ^r= -A + (r)^(r) = V T ^(r), 

where the operator V T depends on r and (x,Xq), and we define the solution operator S 
(6.6) = <W(0) . 



As in (6.5) we can view this as choosing special initial data for tp(0). From now on we only consider such 
data. 

The operator V can be symmetrized 

(6-7) V T :=G- 1 V T G T = {V-V Q ) T -^ M Y, K XV 

j 

with real eigenvalues {A m } and orthonormal eigenvectors {C m } in L 2 (dx dX*), satisfying 

V T C(r) = \ m {T)C{r). 

Therefore V r has the same eigenvalues and the eigenvectors (7™ := G T C, m , which establishes the spectral 
representation 



VIM / m-,T,-)x m )G- 2 dx* 



(6.8) VM;T,-) = y\ m {r) (^(•,r,.),C m )G; 2 dX H ,C m (r) 



We note that the weight G 2 on the co-dimension one surface T 3N 1 appears precisely because the operator 
V is symmetrized by G~ 2 and the weight G~ 2 corresponds to the Eulerian-Lagrangian change of coordinates 
d3T4| 

(6.9) / (i>,( m )G- 2 dX*= [ (i>,C m )dX . 

7 T 3N~1 JJ3N-1 

The existence of the orthonormal set of eigenvectors and real eigenvalues makes the operator V self-adjoint in 
the Lagrangian coordinates and hence the solution operator S becomes unitary in the Lagrangian coordinates. 

6.2. Stability from perturbed Hamiltonians. In this section we derive error estimates of the weight 
functions G when the corresponding Hamiltonian system is perturbed. To derive the stability estimate we 
consider the Hamilton- Jacobi equation 

H(V X 9(X),X) =0 
in an optimal control perspective with the corresponding Hamiltonian system 

X t = V P H{P u X t ) 
P t = -V x H(P u X t ). 

We define the "value" function 

e(x ) = e(x t )- [ h{p si x s )d s , 

Jo 

where the "cost" function defined by 

h(P, X):=P- V P H(P, X) - H(P, X) 
satisfies the Pontryagin principle (related to the Legendre transform) 

(6.10) H(P, X) = sup (P ■ V Q H(Q, X) - h(Q, X)) . 

Q 
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Let 9 be defined by the hitting problem 



0(X O ) = 0(X T ) 



f h(P s ,X s )ds 
Jo 



using the hitting time r on the return surface I 

(6.11) t := inf{i \X Q e I, X t € Ikt > 0} . 

For a perturbed Hamiltonian H and its dynamics (X t ,Pt) we define analogously the value function 9 and 
the cost function h. 

We can think of the difference 9 — 9 as composed by a perturbation of the boundary data (on the return 
surface I) and perturbations of the Hamiltonians. The difference of the value functions due to the perturbed 
Hamiltonian satisfies the stability estimate 

0(X ) - 6(X ) > 0{X f ) - 9(X f ) + f (H — H) h x 9(X t ),X t ) dt 

(6.12) J° 

9(X ) - §(X ) < 9{X T ) - 9{X T ) + J (H- H) (v x fl(l f ),Xt) dt 
with a difference of the Hamiltonians evaluated along the same solution path. This result follows by differ- 



entiating the value function along a path and using the Hamilton- Jacobi equations, see Remark 6.1 and [9]. 
We assume that 



(6.13) 



sup \(H - H)(P,X)\ =0(M~ 1 ), 

(p,x)=(v x e(x t ),x t ),(p,x)=(v x e(x t ) : Xt) 



which is verified in (6.20) for Schrodinger and Born-Oppenheimer Hamiltonians. We choose the hitting set 
as 

(6.14) / := {X e T m \9(X) = 9{X)} 

on which the two phases coincide. Now assume that / forms a codimension one set in T 3Ar and that the 
maximal hitting time r for characteristics starting on / is bounded; the fact that / is a codimension one set 
holds, for instance, locally if \V x{& — &)\ is nonzero. In fact, it is sufficient to assume that there exists a 
function 7 : T 3N -> R, satisfying 7 = C(M _1 ), and such that the set I := {X e T 3N \ 9{X) - 9{X) = j(X)} 
is a codimension one set with bounded hitting times. Then the representation (6.12), for any time t replacing 
r and f , together with the stability of the Hamiltonians (6.13) and the initial data (9 — 9)\i = obtained 



from (6.14) imply that 
(6.15) 



\\e-d\\ L °c =o(m- 1 ), 

provided the maximal hitting time r is bounded, which we assume. 

When the value functions 9 and 9 are smoothly differentiable in X with derivatives bounded uniformly in 



M, the stability estimate (6.12) implies that also the difference of the second derivatives has the bound 

0(M~ 1+S ) , for any 8 > 0. 

|G|- 2 (^,V) with G defined in ( pDT| ). The Born- 
0(M~ 1 ) thus it remains to estimate the weight 



(6.16) \\A X 9-A X 9\\ L , 
Our goal is to analyze the density function p 



1 



Oppenheimer approximation (5.1) yields {ip^) 
function |G|~ 2 . This weight function satisfies the Hamilton- Jacobi equation 

(6.17) H G (Vx log \G\-\X) := V X 9(X) ■ V x log \G\~ 2 + A X 9(X) 



0. 



The stability of Hamilton- Jacobi equations can then be applied to (6.17), as in (6.12 1, using now the hitting 

set 

(6.18) I:={Xe T 3N I log \G{X)\- 2 = log \G{X)\- 2 } 

and the assumption of bounded hitting times r in the hitting problem, and we obtain 

(6.19) || log |G|- 2 - log |Gr 2 || L ~ < C\\H G - H G \\ L - = 0{Nr 1+& ) . 

In this sense we will use that an 0(M~ 1 ) perturbation of the Hamiltonian yields an error estimate of almost 
the same order for the difference of the corresponding densities p — p. 
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The Hamiltonians we use are 



H BO 



2 
2 



Mx),vm(x)) 

+ X {X)-E, 



E, 



based on the cost functions 



ha = E ■ 



\P\ 2 MX),V(XMX)) 
2 



IPI 2 

/ lBO = £; + LL-Ao(X). 

For the Born-Oppenheimer case the electron wave function is the eigenstate $bo- The Born-Oppenheimer 
approximation (5.1), proved in Lemma 6.2 implies that 

(6.20) \\Hs-Hbo\\l- =0(M- 1 ), 



which verifies (6.13|. 



Remark 6.1. This remark derives the stability estimate (6.12). The definitions of the value functions imply 



6{X f )- I h(Pt,X t )dt— y6{X T ) — I h{P u X t )dt 



(i 



(6.21) 



§(X ) S(X ) 

h{P t ,X t ) dt + 6{X f ) - 9(X )+9{X f ) - 9{X f ) 

0(X ) 

h(P u X t )dt+ { d9(X t ) + 9(X f ) - 9(X f ) 



-h(Pt,X t ) + Vx6(X t ) ■ V P H(P t ,Xt) dt + 0(X f ) - 9{X f ) 



< 



(H — H 



<H(VxO(X t ),X t ) 

)(y x e{x t ),x^j dt + 9(x f )-9{x f ), 



where the Pontryagin principle (6.10) yields the inequality and we use the Hamilton-Jacobi equation 

H(V x 9(X t ),X t ) = 0. 

To establish the lower bound we replace 9 along with X t by 9 and X t and repeat the derivation above. 

6.3. The Born-Oppenheimer approximation. The purpose of this section is to present a case when the 
Born-Oppenheimer approximation holds in the sense that \\tp — ^Bo\\L 2 (dx) 1S small. 

We know from Section 6.1.1 that the solution ip t = Stfltpo is bounded in L 2 (dxdX), since iS is unitary 
in the Lagrangian coordinates. This unitary S implies that the integral in the Lagrangian coordinates 
J T 3jv-i (V'ti ipt) dXo is constant in time. We consider the co-dimension one set 



h := {X € 



»3jV 



0>(x),tP(x)} 



•f 3JV-1 



(iP(t,X ),i>(t,X )}dX / 



J3N-1 



dX }. 



where the point values of (ip(X),tp(X)) coincides with its L 2 average. We choose a time t such that X t € 1^, 
and assume that the time r* it takes to hit 1^ the next time is bounded, i.e., 

t* := inf{r | X t G I 4> , r > & X t+T e 1^} = 0(1) . 

We also assume that all functions of X are smooth. 
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Lemma 6.2. Assume that iip = M 1 ! 2 Vtp holds, then there exists initial data for ip such that the L 2 (dx) 
orthogonal decomposition ip = ipo tpQ , where ipo — «$bo for some a G C satisfies 

WKWW L^(dx) 

\\fPo(t)\\mdx) 



= OiM- 1 / 2 ) 



(6.22) 



uniformly in time t, provided the spectral gap condition (5.2) holds, the smoothness estimate (6.29) is satisfied 
and the hitting time r* is bounded. 

Proof. We consider the decomposition ip — ipQ © tpQ, where ipa{r) is an eigenfunction of V{X T ) in L 2 (dx), 
satisfying V(X t )iPq(t) = X (t)iPq(t) for the eigenvalue Xq(t) E R. This ansatz is motivated by the zero 
residual 

(6.23) Kip:=ip + iM 1/2 ViP = 
and the small residual for the eigenfunction 

(n(^o)>o) = o 

M 1/2 V^ = 0{M~ 1/2 ) , 

where 

(6.24) w(X) = (y BO {X),w(X))V B o(X) © Uw(X) 

denotes the orthogonal decomposition in the eigenfunction direction ^bo and its orthogonal complement 
in L 2 (dx). We consider first the linear operator 1Z in (6.23) with a given function Vq and then we use a 
contraction setting to show that Vb = (ip, Vip)/{ip, ip) also works since ||V>o~||L 2 (d:r) is small. The orthogonal 
splitting ip = ipo © tfjQ and the projection II(-) in (6.24) imply 

o = n(^ + 0) 

= n(^ ) + ^o + i M 1 /2(V-K))Vo ± + l nf a x (g-Vo ± )J , 

where the last step follows from the orthogonal splitting 

n((v - v )i£) = (v- v )^ 

together with the second order change in the subspace projection 

ip L (r + At) = n(r + At)(^(t + At)) = H(t)(^(t + At)) + 0(At 2 ) 

which yields H(iPq) = ipQ-, here II(t)- denotes the projection on the orthogonal complement to the eigen- 
vector ip (r). To explain the second order change start with a function v satisfying (ii,$bo(X t )) = and 
f BO (I ff ) = * BO (X r ) + O(At) for g G [t, t + At] to obtain 

n(ff)(n(T + Ar)«-n(r)«) = n(o-)((v,y BO (x T ))v BO (x T ) - (v^ bo (x t+At ))^ bo (x t+At )^ 

= u(o-)o(Ar 2 ) + n(<r) (<«, o(At))* bo (x t )) 

= C(At 2 ) + 0(At)(* bo (X t ) - $> BO {X T ), * BO (X a ))* BO (X a ) 
= 0{At 2 ). 

Let S Tt(J be the solution operator from time a to r for the generator 

'GM -1 / 2 



u^iA/ 1 /2(v-y ) u + in 



2 
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A X (G- 1 «) =: iM^Vw. 



Consequently, the perturbation ip^ can be determined from the projected residual 

^ = -iM^VVo - n(^ ) 

and we have the solution representation 

(6.25) V^(t) =5 T)0 ^(0)- / S T . a IL(llMa)) da . 

Jo 

Integration by parts introduces the factor M -1 / 2 we seek 

/ § Ta Wlijj (a)da= [ iM- 1/2 -^{S T<J )V- 1 Wl^{a)da 
Jo ' Jo da 



iM l/2 lr ( S T,*V~ lun Mv)) da 



(6.26) 



- J lAr 1 / 2 ^— (y- 1 (X a )im$ (<rj) da 

= iM^^V^imMr) - iM- 1/2 5 T! oV _1 n^ (o) 

-J iM-V 2 S Tta ^(v-\X a )WlTp (a)^ da. 
To analyze the integral in the right hand side we will use the fact 

V- 1 = (l + CV-Vo)- 1 [V-(V-V )\y l (V-V )-\ 

which can be verified by multiplying both sides from the left by I + (V — Vb) _1 V — (V — Vq) . A spectral 

decomposition in L 2 (dx), based on the electron eigenpairs {A^, "0fe}^Li and satisfying Vipk — ^k4>k, then 
implies 



v- x n(^ ) = ( z + (v - Vo)- 1 v - (v - v ) ) (v - vb^iWo) 



(6.27) 



^(i + iv-Vo)- 1 [V-(V-Vb)]) (A fc -7 )-Vfe(n(^o),^) 



k=tO 



= ]T(A fe - y )-Vfc(n(^o),^) + oiM- 1 ) 

fc#0 

which applied to the integral in the right hand side of (6.26) shows that H^Hl 2 ^) = 0{M^ 1 / 2 ) on a 
bounded time interval, when the spectral gap condition holds and V'fc are smooth. 

The evolution on longer times requires an additional idea: one can integrate by parts recursively in (6.26) 
to obtain 



S Tt<T IllZipo(a) da 



da da K da 

B := iM~ l l 2 V~ l , U-.^mZMa), 



(7=0 



so that by (6.25) we have 



^(r)=5 T , ^(0)- 



^ d 



^d ,~d 



By choosing 

we get 
(6.28) 



a=0 



(BK{a) - B^(BK)(a) + B^(B^-(BK))(a 
\ da da da 

oc 

^W = -E^o(r), 



)-..) 



a=0 



n=0 
21 



where Bq :— —iM 1 / 2 V x £p and IZo := iM 1 / 2 V lr R. We assume this expansion (6.28) is convergent in 
L 2 (dx) for each r, which follows from the smoothness estimate 



(6.29) 



and ( |6.27 ). 

The next step, verifying that also the non linear problem for Vo works, is based on the contraction obtained 
from 

(i>, (V - A )V> 



Vo-Ao 



0(\\ri\\L*(d*)) 



and that ipQ depends on Vo in (6.25), (6.261 and ( |6.27| ) with a multiplicative factor 0{M 1//2 ) 
Finally, to conclude that | (-0, V 7 ) - 1 = 0{M- 1 ), we use the evolution equation 

±{^) = M- 1 / 2 |G| 2 Im<A|, £) = O(M-i) 



where the last equality uses the obtained bound of ip^ in the first part of (6.22 ). The assumption of a finite 
hitting time r* then implies that | (ip, tp) — 1\ = 0(t* M^ 1 ) = (D(Af~ 1 ), since we may assume that {ip, ip) = 1 
on I^p. □ 

Remark 6.3 (Error estimates for the densities). We have the densities 



(6.30) 
(6.31) 



p s = G s 2 (V,V> 
Pbo = G 



BO 



for the Schrodinger equation, 
for the Born-Oppenheimer dynamics. 



From the stability of the Hamilton- Jacobi equation for log(|G| 2 ) and the estimate \\dxixi{9 ~ 
0(M~ 1+S ) in (|6.16[) we have 



G, 



Gl 



and Lemma 6.2 implies 
(6.32) 

which proves 



BO -rO(M 

(ip,<ip) = 1 + 0(M~ 1 ) 



Ps = Pbo + 0{M 



7. Fourier integral WKB states including caustics 

7.1. A preparatory example with the simplest caustic. As an example of a caustic, we study first the 
simplest example of a fold caustic based on the Airy function A : K — > K which solves 



(7.1) 

The scaled Airy function 

solves the Schrodinger equation 
(7.2) 



■ d xx A(x) + xA(x) = . 
u(x) = CA(M 1/3 x) 



1 

M 



d X xu(x) + xu{x) = , 



for any constant C . In our context an important property of the Airy function is the fact that it is the 
inverse Fourier transform of the function 



Hp) 



7T 



(7.3) 



A(x) = - / e 

7T 



i(xp+p 3 /3) 



dp . 



In the next section, we will consider a general Schrodinger equation and determine a WKB Fourier integral 
corresponding to (7.3) for the Airy function; as an introduction to the general case we show how the derive 
(7.3): by taking the Fourier transform of the ordinary differential equation (7.1) 

(-d xx + x) A(x)e- ixp dx = (p 2 + id p )A(p) , 



(7.4) = 

we obtain an ordinary differential equation for the Fourier transform A(p) with the solution A(p) = Ce lp , for 
any constant C . Then, by differentiation, it is clear that the scaled Airy function u solves ( 7.2 ). Furthermore, 
the stationary phase method, cf. Section [9j shows that to the leading order u is approximated by 

u(i)-c(-ri/ 1/3 j cos(M 1/2 (-x) 3/2 -ir/4) , for x < , 

and u(x) ~ to any order (i.e., 0(M~ K ) for any positive K) when x > 0. The behaviour of the Airy 
function is illustrated in Figure [3J 




Figure 3. The Airy function. 



7.1.1. Molecular dynamics for the Airy function. The eikonal equation corresponding to (7.2) is 

p 2 + x = 

with solutions for x < 0, which leads to the phase 

(7.5) p = 6'{x) = ±(-x) 1/2 , and 9(x) = T \{-xfl 2 . 

o 

We compute the Legendre transform 

9* (p) = xp — 6{x) 

where by (7.5) and —x—p 2 we obtain 



2:t 



We note that this solution is also obtained from the eikonal equation 

P 2 +d p 9*(p) = 0, 

which is solved by 

0*(p) = -p 3 /3. 

Thus we recover the relation for the Legendre transform — xp + 0*(p) = —9(x). 

7.1.2. Observables for the Airy function. The primary object of our analysis is an observable (a functional 
depending on u) rather than the solution u(x) itself. Thus we first compute the observable evaluated on the 
solution obtained from the Airy function. In the following calculation we denote by C a generic constant 
not necessarily the same at each occurrence, 

g(x)\u(x)\ 2 dx = C [ g(x) f f e -M 1/2 ( W/s^M^fW/s) dqdpdx 

JR JR Js. 

= C J J g (M^ a (p - q)) e * Ml/2 (? 3 /3-P 3 /3) dqd p 
= C [ [ g(M 1/2 {p - q)) e ^ 1/2 ((?-p) 3 /i2+( 9 -rt( P + 9 ) 2 /4) dqdp 

q=q-p,p=p+q 

= C J [ g( z M^q)e iMl/2 ^/ 1 ™^ i ) dq^dp 
= C I [ g(t)^ {t3/il2M)+tf2,A) dtdp 



(7.6) 



= C g*A M ( -p 2 )dp 

=d p 0*(p)=x 

= C g* A M (x)\d x p(x) \ dx , 

J — oo 

where 

/ M\ 1 ^ 3 ( I ' M\ 1 ^ \ 
(7.7) A M (x) := ( — j A ( — J a; is the Fourier transform of e -" 3 /(i2M) 

Lemma 7.1. The scaled Airy function Am is an approximate identity in the following sense 
(7-8) ||s * A M - g\\ L 2 (R) < j^\\dlg\\ L 2 {R) . 

Proof. Plancherel's Theorem implies 

M\\g * A M - g\\ L , = M\\gA M - §\\ L * = \\g(e^^ 2M ^ - l)M\\ L , 



1 ■■, ,., 1 

L 2 



<^Hbl 3 slU» = tttII^ 



12" |jr| 12' 
The inequality follows from \e ly — 1| < \y\ which holds for all y G M. □ 



The classical molecular dynamics approximation corresponding to the Schrodinger equation (7.2) is the 
Hamiltonian system 

X=p, i> = -\ 

with a solution X t = —t 2 /4 and the corresponding approximation of the observable 
T Jo T J X t T J_ T 2 /A \p[x)\ 



In this specific case the phase satisfies \p{x)\ — Ix] 1 / 2 and \d x p\ — \x\ 1 / 2 /2 , and hence the non-normalized 
density is in this case equal to 2|dt E _p|. Equation (7.6) and Lemma 7.1 imply 

| / g\u\ 2 dx- f gd x p(x) dx\ = 0(M~ 1 ) 

JM JUL 

and consequently for two different observables g\ and g 2 we have that Schrodinger observables are approxi- 
mated by the classical observables with the error 0(M~ 1 ) 



(7.9) 



J R 9iH 2 dx _ J R gi\d xx 9\dx 
J R g 2 \u\ 2 dx J R g 2 \d xx 9\dx 



OiAf- 1 ) 



using d x p(x) — d xx 9(x). The reason we compare two different observables with a compact support is that 
J m u 2 (x) dx = 00 in the case of the Airy function. 
We note that in (7.6) we used 



1 



- p*) = 9* (p) - 9* (q) = (p - q)d p 9* l-(p + q)\+-9 



1 



1 



(p + ?) ) ( 1(p-q) 



which in the next section is generalized to other caustics. For the Airy function there holds 



1 



^9* [-(p + q) =-- 



7.2. A general Fourier integral ansatz. In order to treat a more general case with a caustic of the 
dimension d we use the Fourier integral ansatz 



(7.10) 

and we write 



4>(X, x)e 



-iAI 



2 e ( x,x,p) d p 



X=(X,X), P=(P,P) 

d N 

A /' ^\\'/". X-P= X j P j 

3=1 J=d+1 

Q(X,X,P) = X ■P-0*{X,P) , 



based on the Legendre transform 



x 



d*{X,P) =mm(X-P-9(X,X) 



If the function 8*(X, P) is not defined for all P € M. d , but only for P e U C M. d we replace the integral over 
~R d by integration over U using a smooth cut-off function x(-P). The cut-off function is zero outside U and 
equal to one in a large part of the interior of U, see Section [7.2.3 The ansatz (7.10) is inspired by Maslov's 
work 25], although it is not the same since our amplitude function 4> depends on [X,X,x) but not on P. 
We emphasize that our modification consisting in having an amplitude function that is not dependent on P 
is essential in the construction of the solution and for determining the accuracy of observables based on this 
solution. 



7.2.1. Making the ansatz for a Schrodinger solution. In this section we construct a solution to the Schrodinger 
equation from the ansatz (7.10). The constructed solution will be an actual solution and not only an 
asymptotic solution as in |25j. We consider first the case when the integration is over M. d and then conclude 
in the end that the cut-off function x(P) can be included in all integrals without changing the property of 
the Fourier integral ansatz being a solution in the A-domain where X — V p9*(X, P) for some P satisfying 

x(P) = 1- 

25 



The requirement to be a solution means that there should hold 

(7.11) 

= (H- E)& 

= J Q|V*0*(X, P)\ 2 + l\P\ 2 + V (X) E^J <j>(X, x)e -iM^e ( x,x,P) dP 
- j (iM-^Vtt ■ V A 9* - V t <j> ■ P + 1^9*) - (V - W + ^Az^j e-u^et*,*,*) dA 

Comparing this expression to the previously discussed case of a single WKB-mode we see that the zero order 
term is now A^(9* instead of Ax9 and that we have —V x 4> • P instead of • However, the main 

difference is that the first integral is not zero (only the leading order term of its stationary phase expansion is 
zero, cf. (9.1 )). Therefore, the first integral contributes to the second integral. The goal is now to determine 
a function F(X,X,P) satisfying 



1 



IV *0 



* |2 



1 



P\ A +Vo(X)-E)e 



(7.12) 



and verify that it is bounded. 

Lemma 7.2. There holds F = Fq + Fx where 

Fo = IJ2 9 x*xi V (V P 9* (P)) d PiPi 9* (P) 



-iM 1/2 e{x,x,P) 



iM~ x l 2 I F(X,X,P)e- iMl/2@ (*>*> p UP, 



dP 



Fi = iM- 1/2 



i f i 



o Jo 



*(1 - t)dp« [dx'xix*Vo (Vp6»*(P) + stS9*(P)) dp.p.V p9*(P) 



dt ds . 



Proof. The function 9*(X,P) is defined as a solution to the Hamilton- Jacobi (eikonal) equation 
(7.13) ^\V jt 6*(X,P)\ 2 + ^\P\ 2 + Vo(x,Vp6*(X,P))-E = 

for all (X,P). Consequently, the integral on the left hand side of ( 7.12| ) is 

(V (X, X) - V (X, V P 9*(X, P)) e -iMV*{x.p-6* { x,P)) d p 

Let Po(X) be any solution to the stationary phase equation X — V p9*(X, Po) and introduce the notation 

Q'(X,X,P) :=Vp9*(X,Po)-P-6*(X,P). 

Then by writing a difference as V(yi) — V{y2) — Jq d y V(i/2 + — V2))dt ■ (yi — y 2 ), identifying a derivative 
dp and integrating by parts the integral can be written 

/ (v o (X,Vp9*(X 1 P Q ))-Vo(X,Vp9*(X,P))e- lMl/20 '^^^dP 
= 11 Y, d x^o{Vpe*{P)+t [Vp9*(P )-Vp9*(P)])x 

JO JR d t 

x (dp,9*(P ) - dp t 9*(P)) e -^B'(x,x,P) d p dt 
= -iM~ 1 ' 2 



= iM- 



J2 d x* v o {Vp9*(P) + t [Vp9*(P ) - Vp9*(P)]) dhe -iM^e'(x,A,P) d p dt 

i 

1/2 t I Y J dp.d jti V Q {Vp9*{P)+t [Vp9*(P )-Vp9*(P)])e- lMl/2e '^^> p UPdt. 

JO JTR d i ' 



2(> 



Therefore the leading order term in F =: Fq + F\ is 



F Q := I Y j {l-t)d xixi V (Vp9*{P))dp i p i 9*{P)dt 



Denoting 89* (P) — Vp9*(Po) — Wp9*(P) the remainder becomes 



iM~ 1 ' 2 



I f E i B *'**Vo (Vpr(p)) -a M vb (v P e*(p) + ****(£))] 



(1 - t)dpj p%9*(P) e -^e'(x,x.p) d p dt 



= iM- 1 ' 2 



I I I T, t ( 1 -^ d x^x^o(^pO*(P) + stS9*(P))d P3P ,9*(P) 
Jo Jo JR d iJ:k 

< (dp k 9*(P ) - d pk 9*(P)) e -iM^e'(x,x,P) d p dtds 

h [ [ I E^ 1 "*)^* [dx*xix*V (Vp9*(P) + stS9*(P))dp j p i 9*(P)] 
V1 Jo Jo Jm d TTu 



xe -iM^e'(x,x,P) dPdtdSi 
hence the function F\ is purely imaginary and small 



F 1 = iM~ 1 / 2 / / ^t{l-t)dp h d^xj^Vo (Vp9*(P) + st 89* (P)) dp. P yp9*(P) 



dt ds . 



and 
(7.14) 



2ReF = E 5 *'*^ (VpO*(P)) d P] p,9*(P) 



a 



The eikonal equation (7.13) and the requirement that (H — E)<& = in (7.11) then imply that 
= 



(7.15) 



iM~ 1/2 [ V x (f>-V x t 



1 



- V x <f> ■ P + -</, (A x 9* - 2F(X, P)) 
e -iM 1/2 e(x,x,P) d p 



The Hamilton-Jacobi eikonal equation (7.13), in the primal variable (X,P) with the corresponding dual 
variable P, X), can be solved by the characteristics 



(7.16) 



using the definition 



The characteristics give 



X = P 

P = -V X V (X,X) 
X = -P 

P = V X V (X,X), 

V X 9*{X,P)=P 
Vp9*{X,P) = X. 



-<l> = V x 4,.V x 9*-\7 x ct>-P, 
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so that the Schrodinger transport equation becomes, as in (3.12), 



(7.17) 



and for tp = G(j) 
(7.18) 



with the complex valued weight function G defined by 

d 



(7.19) 



dt 



log G t 



A x 9*(X t ,P t )-F(X t ,P t ). 



This transport equation is of the same form as the transport equation for a single WKB-mode, with a 
modification of the weight function G. 



Differentiation of the second equation in the Hamiltonian system (7.161 implies that the first variation 
dP t /dX satisfies 

A / f)pi\ r)P k 



dt \dX c 



dXn 



which by the Liouville formula (3.14) and the equality 



2ReF = J2dxixiVodp j p i 9* = Tr(^\v Y v V dp r J>\ 



in (7.14) yields the relation, 



(7.20) 



-2 / * Re_Fdt' _ Q 



det 



dP t 



dX 



for the constant C 



j det ^prl- We use relation (7.20) to study the density in the next section. 



Remark 7.3. The conclusion in this section holds also when all integrals over dP in M. d are replaced by 
integrals with the measure x(-P) dP. Then there holds 2ReF = . dxixjVdpi{xdpj6*). We use that the 
observable g is zero when the cut-off function \j is n °t one, see Section 
how to construct a global solution by connecting the Fourier integral so 



7.2.3 In Section 7.2.5 



we show 
rborhood 



lutions, valid in a neig 

where detd(X)/d(P) vanishes (and x(P) = 1)> to a sum of WKB-modes, valid in neighborhoods where 
detd(P)/d(X) vanishes (and x(P) < !)• 

7.2.2. The Schrodinger density for caustics. In this section we study the density generated by the solution 



*(X,x)= [ <t>(X,x)e- iMl/2 ( jt - p - e "(*' p »dP. 



The analysis of the density generalizes the calculations for the Airy function in Section 7.1.2 We have, using 
the notation g for the Fourier transform of g with respect to the X variable, and by introducing the notation 
R= \{P + Q) and S = P-Q 
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l{X)\*(x,X)\*dxdX = I g{X){^ct>) e iM^{x.P-e-{x,p)) e -iMVHX-Q-e%&,Q)) d p d Q d x 



■MX) 



g(X, M X ' 2 S) e iM^(0*(x,Q)-e' (x,f)) d p d Q d x 



= I g(x,M 1/2 5)e lMl/2 B (5 ' v - fi)3e * ( ^' A+75/2) : 



(7.21) 



x e 



iM 1/2 S-Vp0*(X,R) 
d/2 



dSdR dX 



—J M~ 1/2 J g*A M {X,Vp6*(X,R))dRdX 



Q^j ' M-V3 J~g*A M (X,X) 



=x 



det 



a(p) 



dX". 



In the convolution g * Am, the function Am, analogous to (7.7), is the Fourier transform of 

e »i( U -V P ) 3 9'(X,P) 

P=R+ 7 UJ 

with respect to w € M'' and the integration in X is with respect to the range of V p9*(X, •). As a next step 
we evaluate the Fourier transform and its derivatives at zero and obtain 



A M (X)dX = l, / X l A M (X)dX = 0, 

X l X 3 A M (X) dX = , M f X l X J X k A M (X) dX = 0(1). 



Here we use that both differentiation with respect to (w ■ V P ) 3 and 9*(X, R + jut) yield factors of u> which 
vanish. The vanishing moments of Am imply that 



(7.22) 



\\g*A M -g\\ L2{dJt) =0(M l ) 



as in (7.8), so that up to 0{M - 1 ) error the convolution with Am can be neglected. 



7.2.3. Integration over a compact set in P. In the case when the integration is over U C M. d instead of R rf , 
we use a smooth cut-off function x(-P), which is zero outside Li and restrict our analysis to the case when 
the smooth observable mapping P i-> g(X, Wp9*(X, P)) is compactly supported in the domain where x is 
one. In this way g(X, V p9* (X, P)) is zero when V P %(P) is non zero. The integrand is thus equal to 

(g(X) (<M))x(^MQ) 

and we use the convergent Taylor expansion 

f, A k=0 



Then the observable becomes 



oo . 

{2^)- d / 2 M- 1 ' 2 Y, / (a k (M- i A jt ) k g)*A M (x,Vpe*(X,R^ dRdX . 



2!) 



As in ( 7.22[ ) we can remove the convolution with Am by introducing an error 0(M 1 ) and since for k > 
we have ak(R)g{X , V p9* (X , R)) — and ao = 1, we obtain the same observable as before 



oo 

E 

k=0 



a k (M^Axfg) * A M (X, V p 6* (X,R)) dRdX 
jP f (a k (M~ 1 A x ) k g) (x, V P 9* (X, i?)) dRdX +0(M- 1 ) 

J g(x,Vp9*(X,R)) dRdX + C(M _1 ) 
dX +0(M- 1 ). 



~g{x,x) 



dot 



g(^) 



7.2.4. Comparing the Schrddinger and molecular dynamics densities. We compare the Schrodinger density 
to the molecular dynamics density generated by the continuity equation 

= div(pV6») = Vp ■ V6> + pdiv(V6>) = p + pdiv(V6») 

which yields the density 

e -/div(ve)dt _ 

We have P = V(9, so that §jj^j = dxx9- The Liouville formula (3.14| implies the molecular dynamics 
density 

dX bo 



(7.23) 



Pbo 



= e -I^ iv (ve)dt' =det 



<9X 



t,BO 



The observable for the Schrodinger equation has, by (7.21), the density 



dot 



9(P) 



We want to compare it with the molecular dynamics density pbo- The convolution with Am gives an error 
term of the order (D(M~~ 1 ) : as in (7.8), and following the proof of Theorem 5.1 for a single WKB-state in 



Section [6] (now based on the Hamilton- Jacobi equation (7.13), the Schrodinger transport equation (7.17) 



Born-Oppenheimer approximation Lemma 



and the definition of the weight G in (7.19)), the amplitude function satisfies, by (7.18) and (7.19) and the 

e /2Rc F-A x e* dt 



6^ 



so that by (7.20) 



(g((/>, (f>)) * A M |det 



■ = |G| 2 (V,V> 
d(P) 



d(xy 



(g{<P,cj>})\ det 



d(P) 



(7.24) 



d(xy 

ffe /2RcF-A^- dt | det 

d(X Q ) 



= g I det 
= 5 I det 
= 9 | det 



0(P) 

g(jg) | 

0(*o) | 
d(X) 1 



I det 



0(M~ 

d(P) 
d{X) 

d(X ) 



d(X) 
d(X ) 
d(X) 1 

-0{M- V ) 



det 



0(M- r 



det 



d(X)< 



OiM- 1 ). 



When we restrict the domain to Li with the cut-off function \ as i n Remark 7.3 we use the fact that 
g(X, Wpd*(X, P)) is zero when V px(P) 1S non zero an d obtain the same. The representations (7.24) and 
(7.23) show that the density generated in the caustic case with a Fourier integral also takes the same form, 
to the leading order, as the molecular dynamics density and the remaining discrepancy is only due to 9* = 6g 
and 9* = # BO being different. This difference is, as in the single mode WKB expansion, of size 0(M~ 1 ) 
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which is estimated by the difference in Hamiltonians of the Schrodinger and molecular dynamics eikonal 



equations. The estimate of the difference of the phase functions uses the Hamilton- Jacobi equation (7.13) 



for Og(X,P) and a similar Hamilton-Jacobi equation for 9^ (X,P) with Vb = Abo + 0{M 1 ) replaced by 
Abo- The difference in the weight functions log(|G(X, P)| -2 ) is estimated by the Hamilton-Jacobi equation 



{w x 9*(X,P)-W x -W x V (X,X)-W P ) log\G s (X,P)\- 2 -A A e*(X,P)+ReF(X,P) 



0, 



where ReF is given in (7.14), and by the similar Hamilton-Jacobi equation with Vq = Abo + 0{M 1 ) 



replaced by Abo and 0$ by # B o- 

7.2.5. A global construction coupling caustics with single WKB-modes. We use a Hamiltonian system to 
construct solutions to the Schrodinger equation. Given a set of initial points Xq £ K 3Ar the solution paths 
{(X t ,P t )e R 6N | < t < oo, H(P , X Q ) = E} of the Hamiltonian system 

X t = VpH{P u X t ) 

P t = -V x H{P u X t ) 

with a smooth and bounded Hamiltonian H(P,X) generate a 3iV-dimensional manifold called Lagrangian 
manifold. The Lagrangian manifold defined by the tube of trajectories is defined by the phase function 
0(X) that plays the role of a generating function of the Lagrangian manifold. Thus we seek a function 
9 : U C M. 3N — > K such that P t = V xd{X t ). We show that there exists a potential function 9 by determining 
an equation that preserves the symmetry for the matrix Q t , defined as Q^{X) := dxiP l (X) and Q\ J := 
Q ij (X t ). The relations P* = P l {X t ) and Q\ 3 := d X3 P l (X t ) imply 



PI = jP l (X t ) Y^XU^ PI =J2x{Qf =J2d P3 H{P(X t ),X t )Q 



t ! 



so that 



> x ^J29piH(P(X t ),X t )Q 



+ Y / d P] piH(P(X t ),X t )dx^P l Q l 



=T. 3 X?d xkxJ Pl=Y. 3 Xid x3xk P?=Ql k 



and 



dx^Pl = -d xh (d X iH(P(X t ),X t )) = -d X i X HH(P{X t ),X t )) -^d X i Pt H(P{Xt),X t ) d x *P J 



together with the symmetry of Qt show that 

Qf = -d X i X uH{P t ,X t ) - dpi P iH(P u X t )QfQ\ j 

(7.25) 



J2dp 3Xk H(P t ,X t )Q? -J2d Pixi H(P t ,X t )Q 



Since the Hamiltonian is assumed to be smooth it follows that the right hand side in (7.25) is symmetric and 



3.1.3) 



oo and such points 



thus the matrix Q t remains symmetric if it is initially symmetric. Hence there exists a potential function 
9{X) such that P(X) — Wx9(X) in simple connected domains where Q is smooth. The function Q may 
become unbounded due to the term dpjpiH Q kl Q l; ' , even though H has b ounded third derivatives. Points 
X t at which |Tr (Qt)\ — oo satisfy, by Liouville's theorem (see Section 
are called caustic points. 

The same construction of a potential works for the local chart expressed as X = X(P) instead of P = 
P(X). In fact any new variable X (not including both X 1 and P % for any i), based on 37V of the 6iV 
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detff 



variables (X, P), and the remaining variables 37V variables, P, represent the same Hamiltonian system with 
the Hamiltonian H(P,X) := H(P,X). The Lagrangian manifold is defined by P — V X 8(X) in the local 
chart of P-coordinates with the generating (potential) function 9(X) defined in domains excluding caustics, 



i.e., where det 



dX„ 

ax t 



< oo. Maslov, 



, realized that a Lagrangian manifold can be partitioned, by changing 

coordinates in the neighborhood of a caustic, into domains where P = \7 X 8(X) is smooth. He used the 
generating (potential) functions 8 to construct asymptotic WKB solutions of Fourier integral type. A sketch 
of this general situation is depicted in Figure [4] In previous sections we have described global construction 
of solutions in a simpler case without caustics, i.e., P t = X?x8(X t ) holds everywhere. In this section we 
describe the global construction of WKB solutions in the general case when caustics are present. 





Figure 4. The left figure depicts a graph of the molecular dynamics potential X(X) in 
the case which exhibits caustics at X = a and X — b for a given energy E. The right 
figure shows a general case of the Lagrangian manifold with two caustic points X = a and 
b and its covering with charts Ui. In the charts U, i = 2,4 the manifold is defined by 
P = Wx8i(X) and the solution to Schrodinger equation is constructed by simple WKB 
modes. The caustics belong to the charts Ui, i = 1, 3 and in this case the manifold is defined 
by X = S7x8i(P) and the solutions are given by the Fourier integrals. 



We see that the weight function G, in (3.14), based on a single WKB-modc (3.1) blows up at caustics, 
where det(d(X)/d(P)) — 0, and that the weight function G in (7.17) for the Fourier integral (7.10) blows 



up at points where det(d(P)/d(X)) vanishes. Therefore, in neighborhoods around caustic points we need 
to use the representation 8*(X,P) of the phase based on the Fourier integrals, while around points where 
det(d(P)/d(X)) vanishes we apply the representation 8(X,X) based on the Legendre transform, as pointed 
out by Maslov in J2S] and described in the simplifying setting of the harmonic oscillator in [TT] . 

One way to make a global construction of a WKB solution, which is slightly different than in 25 , is to 
use the characteristics and a partition of the phase-space as follows, also explained constructively by the 



numerical algorithm 8.2 in the next section. Start with a Fourier integral representation in a neighborhood 
U of a caustic point, which gives a representation of the Schrodinger solution $ in U. Then we use the 
stationary phase expansion, see Section [9J to find an asymptotic approximation $ (accurate to any order 
TV € N) at the boundary points X of U as a sum of single WKB-modes with phase functions 8j 



x(P)e 



-iAI 1/2 (X-P- 



- lMl " e ^ x U {X)+O{M- n ) 



where each phase function 8j(X) := X ■ Px,j — 9* (X , Px.j) corresponds to a branch of the boundary and 
the index j corresponds to different solutions Pxj of the stationary phase equation X = Wp8*(X,Px)- 
The single WKB-modes 4>(x,X)e tMl,2e ( x ^ are then constructed along the characteristics to be Schrodinger 
soluti ons in a domain around the point where det(d(P) /d(X)) vanishes, following the construction in The- 
orem 3.1 using the initial data of $ at dU. We note that the tiny error of size 0(M~ N ) that we make in 
the initial data for <f) also yields a tiny perturbation error in <j> of size (D(M~ N ) along the path, due to the 
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assumption of the finite hitting times. A small error we make in the expansion therefore leads to a negligible 
error in the Schrodinger solution and the corresponding density. 

When a characteristic leaves the domain and enters another region around a caustic we again use the 
stationary phase method at the boundary to give initial data for (X, P,4>,G). When the characteristic 
finally returns to the first boundary dhi, there is a compatibility condition to have a global solution, by 
having the incoming final phase equal to the initial phase function in C . We can think of this as trying to 
find a co-dimension one surface I in M? N where the incoming and outgoing phases are equal. First to have 
one point where they agree is possible if we restrict the possible solutions to a discrete set of energies E, i.e., 
the eigenvalues, and therefore the compatibility condition is called a quantization condition. Then, having 
one point where the difference of the two phase function is zero, we can combine this with the assumption 
that the Lagrangian manifold generated by the characteristics path (X t ,Pt) is continuous: the two phases 

have the same gradient on J, since (X,P) = (X,V x 0{X)) = ((X, Vp6*{X, P)), (V X 9*(X, P), P)) so the 

phases are C 1 . In this way we define the (X, P,<j>,G) globally, for the eigenvalue energies E. To evaluate 
observables we use a partition of unity to restrict the observable to a domain with a single representation, 
cither a Fourier integral representation for a caustic or a single WKB-mode when det(d(P)/d(X)) = 0. 

8. Numerical examples 

In order to demonstrate the presented theory we consider two different low dimensional Schrodinger 
problems. For both of these problems we show that there exists a Schrodinger eigenfunction density which 
converges weakly to the corresponding molecular dynamics density as M — > oo with a convergence rate 
within the upper bound predicted in the theoretical part of this paper. 

8.1. Example 1: A single WKB state. The first problem we consider is the time-independent Schrodinger 
equation 

(8.1) H$ := (--^dxx +V\ * = E9 

with heavy coordinate X € (— 7r, n] and two-state light coordinate x € {x-, x+ }. Periodicity is assumed over 
the heavy coordinate, &(X, x) — <!>(X + 2tt,x), and the potential operator V is defined by the matrix 

V(X) \V{X)e(X) + c 

\V{X)e{X)+c 

where we have chosen V(X) = — 2cos(A) + cos(4A), e(X) = 1 + X 2 and c to be a non-negative constant 
relating to the size of the spectral gap of V. The action V<I> is thus defined by 

*(X,x+) 

For each X the potential matrix ( |8.2[ ) gives rise to the eigenvalue problem 

V(X)v = X±(X)v 

with the eigenvalues 



(8.2) V(X) 



(V$)(X,-) = V(X) 



X ± (X) - l - (v{X) ± Sgn(X)^V(Xy+4(V(X)e(X)/2 



of 



where Sgn(A) = ±1 as defined below. When constructing the molecular dynamics density for this problem 

PmdP0= ^-a ( x))' 

one has to determine on which of the two eigenfunctions X± to base this density. When c = the difficulty 
that the eigenvalue functions A + and A_ can cross is added to the problem. In order to determine the 
continuation of eigenvalue functions at the crossings we introduce a function Sgn(A) which is a sign function 
with Sgn(— 7r) = 1 that changes sign at points where 

y(l) 2 + 4^y(X)e(X) + C j =0. 
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Since this situation can only occur when c = 0, it is possible to set 



Sgn(A) := S gn(V(-n))sgn(V(X)) . 



See Figure [5] for a typical eigenvalue function crossing, which makes the function A± : K — > K smooth (in 
contrast to the choice Sgn =1). 

To solve ( |8.1[ ) numerically, we use the finite difference method to discretise the operator M on a grid 
{Xj}f =1 x {x-,x + } with the step-size h = 2-k/N and Xj — jh. The discrete eigenvalue problem 



is solved for the 10 eigenvalues being closest to the fixed energy E and a molecular dynamics approximation 
of the eigensolution is constructed by 

$mb(X,x) := ^p M B(X)e lMl/20 ^v(X 1 x), 
where v(X, •) is one of the eigenvectors of V(X) and 



(8.3) 



O(X) := / ^/2(E 1 - A(s)) ds 
Jo 



is approximated by a trapezoidal quadrature yielding . Thereafter a Schrodinger eigensolution which 
is close to the molecular dynamics eigensolution is obtained by projecting <&md onto the subspace spanned 
by {T}/ =1 as described in Algorithm^] By denoting p$ w (X) = ($( h \$( h )) and pmd(X) = ($md,*md), 
the observables Qi{X) = X 2 and g2(Xj = V(X) are used to compute the convergence rate of 



(8.4) 



X! w gi (x) PMB (x) dx - /* w 9i(x) P9W (x) dx 



j: Tr g l (X) PMU (X)dX 



as M increases. Further details of the numerical solution idea are described in Algorithm [T] 

Plots of the results for the test case with the spectral gap c — 5 and E — 0, and for the test case with 
crossing eigenvalue functions when c = and E = 1.2 are given below. Most noteworthy is Figure [8] which 
demonstrates that the obtained convergence rate for (8.4) is 0(M~ 1 ) for both scenarios. 



Plot of eignvalue functions 



10 


i 
i 




1 


i 
i 




1 

1 


8 
6 


i 

V' \ 

\ 


X (X) 
S(X) 


f 

!'\! 

i 


4 

2 


V 

\ 

\ 

\ 


/ 











-2 








-4 
-6 


/ ~ 

/ 




\ 

V 



-3-2-10123 
X 



Plot of eignvalue functions 



5 









































X (X) 










/•> 

> / » 
i / i 




S(X) 




'\ ' 
' \ ' 






x / > 
~ \ 








\ ' 
' \/ 




\ 


/ 

1 




\ 






i 






\ 

I 


i 




_ / 


1 




1 


\ 




# / 
i v ' 


t 


A 


/ 
1 


\ * * 












( 

/ 












/ 

J 






-3 


-2 


-1 





1 


2 


3 



Figure 5. Left plot: Eigenvalue functions when c = 5. There is a spectral gap which 
makes the sign function constant S = 1. Right plot: Eigenvalue functions when c = 0. The 
eigenvalue functions exhibit crossing, consequently the function S changes its sign from ±1 
to =pl at the crossing points. 
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Algorithm 1 Algorithm for problems in Example 1 



Input: Energy E; potential functions V, e and c; mass M; number of grid points N and grid {Xi}fL 1 . 
Output: Schrodinger projection density p$(h). 

1. Construct the discrete operator 'MM 1 ' from (8.1 1 using finite differences and solve the eigenvalue problem 
for the 10 eigenvalues being closest to E by using MATLAB eigs(H,10,E). 

2. Sort the eigenvalues and eigenvectors by distance from E and keep only the EiS which are less than 
M~ x / 2 away from E. Let J be the number of kept eigenvalues and Eq the eigenvalue closest to E. 

3. 

for i — 1 to N do 

Solve the eigenvalue problem 

V(X i ,-)v±{X i ,-)=\±(X i )v±(X i ,-), 



where V is the matrix defined in (8.2). 
end for 

4. Construct the molecular dynamics density according to the formula 

(E - X(X))- 1/2 



Pmd(X) = 



/pyhr] (Eo-HX)) 



-1/2 



dX 



where we choose X(X) above from the two eigenvalues X±(X) by the criterion that the eigenvalue chosen 
must fulfil || Alloc < Eq. 

5. Construct a discrete molecular dynamics approximation to the eigenfunction 

(8.5) $md (X, x) = y/pMB(Xy Ml/2e ^v(X, x) , 
where v{X,x) is one of the eigenvectors v±, 

(8.6) e(X):=/ y/2(Ex-\(a))da, 



and we approximate by a trapezoidal quadrature 0^. 

6. Project the molecular dynamics solution $md onto the eigenspace {T. t }/ =1 , J < 10 by Algorithm[| to 
obtain a projection solution $W. 

7. Derive the Schrodinger projection density by 
for i = 1 to N do 

p 9W (Xi) = \^ h \X ijX _)\ 2 + \& h XXi,x+)\* , 

end for 

and scaling p^ h ) = p^w /\\p^w ||. 



8.2. Example 2: A caustic state. Next, we consider the one dimensional, time independent, periodic 
Schrodinger equation 

1 



(8.7) 



2M 



dxx+V)<S> = E<$>, X G (—2"/E, 2y/E) 



with V(X) = X and E — 1. The cikonal equation corresponding to (8.7) is 



1 



-P z 



V{X) = E. 
5quation to 

oo to the density generated from a solution of (|8.7| 



As in Example 1, we would like to use the eikonal equation to construct a numerical approximate solution 
of (8.7) whose density converges weakly as M 
The molecular dynamics density corresponding to this eikonal equation becomes by (3.19) peo = C{E 
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Algorithm 2 Projection algorithm 



Input: Mass M; wave solution $; eigenvalues {Ei}{ =1 and corresponding eigenvectors {Tj}^. 
Output: Schrodinger projection wave solution <&( h \ 

1. Organize eigenvalues by multiplicity by a numerical approximation. Construct a J x J, zero matrix A 
which keeps track of multiplicity relations as follows: 

for i — 1 to J do 
for j = i to J do 

if \Ei - Ej\ < M" 3 / 4 then 

Consider eigenvalues equal since the expected spectral gap is 0(M~ 1 ^ 2 ), and store this relation by 
if Akj — for all k < i then 

Set Aij = 1. 
end if 
end if 
end for 
end for 

2. For vectors b £ {0, 1} J , define the projection 

j,k=l 

and, letting p and p^h.b) denote the densities generated by $ and $>( h > b ' respectively, set 

b* = arg min \\p - p^h.b) ||. 
be{o,i} J 

3. Return the projection $W := $( fc » 6 *). 




Figure 6. Plot of the MD density pmd and the Schrodinger projection density p^m in the 
case c = 5 and E = for the two different masses M = 90 (left plot) and M — 724 (right 
plot) illustrating the convergence of the densities. 



V{X))~ X / 2 . The density peo g° es to infinity at the caustics X — V^ 1 (E) = ±y/E and the approach in 
Example 1 does not work directly. We will instead construct the numerical approximate solution using the 
stationary phase method as outlined below based on the WKB Fourier integral ansatz. 
By the Legendre transform 

6*(P) = min (XP-6(X)) 
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Density comparison when M=724 



Density comparison when M=5792 
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Figure 7. Plot of the MD density Pmd and Schrodinger projection density p$(h) in the 
case c = and E — 1.2 for the two different masses M = 724 (left plot) and M = 5792 
(right plot). 
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Figure 8. Left plot: Plot of the observable density errors given in (8.4) with an eigenvalue 
gap, when c = 5 and E = 0. Right plot: Plot of the observable density errors given in (8.4 1 
with an eigenvalue crossing, when c = and E = 1.2. 



an invertible mapping between the momentum and position coordinates fulfilling X = Vp9*(P) is con- 
structed. Using equation (8.8), one sees that V p6*(P) = V^ 1 (E — P 2 /2). Since 6*(0) = 0, one can derive 
that for this particular choice of V 



r p 

9*(P) = / ^E-s 2 /2ds = 
Jo 



E 
72 



sin 



P 

7TB 



/2E V 2E 



In neighbourhoods of the caustics [— 2E 1 / 2 , — Xq) and (Xrj, 2E 1 / 2 ], we construct the approximate solution 

by 

u(X) 



where u is the inverse Fourier transform 



$(X) = 



y/WxV{X)\ 



u(X) := I e ^>\-xp+e* { P)) dp 

-2V~E 
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and Xo E (—V 1 {E), V 1 (E)) is a value yet to be chosen. In the region (— Xq, Xq) the approximate solution 
is constructed by 



(8.9) 

Here 
(8.10) 



*P0 = c 



u(X) 



(E-ViX)) 1 /* ' 



with, according to the Legendre transform, 6(X) := X y / 2(E - V(X)) - 9* \ x/2{E - V(X))J and ij)± de- 
termined by the stationary phase method: 

1. Set P(p) = P a +p with P Q = y/2(E - V(X Q )) and let 



Y{p) :=sgn(p)W2- 



-x(p a +p) + e*(p Q + P ) + 9(x Q ) 



d PP 6*(Po) 



using 



(8.11) 



9(X) := Xy/2(E - V(X)) - 9* ^2(E - V(Xjf) 



and determine its inverse p{Y) in a neighbourhood of Y = by computing (pi,Y(j?i)) on a grid 
around p — and, for fc > 3, fit a 3fc + 1th degree polynomial to the values (Y(pi),pi) using the 
method of least squares. 
2. Evaluate the stationary phase expansion 



u(X ) 



E 



o i7rsgn(d P p0*(\Po))/4 



3.12) 



■p =± v / 2(E-V(X )) 



dp P 9*(P ) 



-1/2 



-iM 1/2 8(X Q ) 



X E^T^ f*f^i'^ 1 <A > i I'M 



3=0 



+ 0(M- l/2 ) 



to obtain 



t (X ") = e lMl/29 ^)(V+ + C(M- fe / 2 )) + e- lMl/2f, ( x °)(^_ + 0(M- fc / 2 )) , 



where 



V>± := 



ol7 rsgn(a P p9*(±Po))/4 



-d PP 6*(±p ) 



fe=0 



A;! 



<9- 



yy 



y=o 



The constant C in ( |8.9[ ) is chosen so that the wave solution parts are continuous at the gluing point, 
&(±Xq) — $(±Xq~). It is most easy to determine C when Xq is chosen so that |u(Xo)| is at a local 
maximum; see Figure [9] for an illustration of the gluing procedure. 

At the end a Schrodinger eigenfunction solution $W is obtained by projecting $ onto the space spanned 
by a set of eigensolutions to the discretized version of the Schrodinger problem, {Tj}^ =1 , as is described in 
Algorithm [2j 

Two convergence results are needed to make the method work. First, the density generated from the 
stationary phase based on the approxmiate solution p(X) :— |<I > | 2 (X)/||<I , || 2 , must converge weakly to the 

for an 
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Schrodinger projection based density p^( h )(X) :— |$(' l )| 2 (X)/||<£>(' l )|j 2 , as M — > oo; see Figure 
illustration of how these functions converge. Second, p$(f») must converge to the molecular dynamics density 
Pmd(-^0 := C(E — V(X))^ 1 / 2 as M increases; see Figure 
A numerical test of the convergence rate of 
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(8.13) 



J-^9i(X)pmd(X) dX f^- o 9i(X)p mW (X) dX 



J11%92(X) P mb(X) dX f^g 2 (X)p^ h) (X) dX 
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— W 

— u(X) 
o gjuirgpoirts 




FIGURE 9. Plot illustrating the gluing procedure of the functions u(X) and u(X) at the 
points ±Xq. 
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Figure 10. Comparison of the approximate solution based density p and the Schrodinger 
projection based solution p^h) for M = 200 (left plot) and M = 800 (right plot). 



as M increases is illustrated in Figure 12 for the observables 



(l.5-Xf(1.5 + Xni + e- x2 ) (1.5-X) 6 (1.5 + X) 6 (1-X 2 + X 4 ) 
(8.14) g^X) = j--gi2 and 92{X) = j-gia ■ 

Further details of the solution procedure in Exampe 2 are given in Algorithm [3] 

9. The stationary phase expansion 

Consider the phase function X ■ P — 6*(X,P) and let Po(X) be any solution to the stationary phase 
equation X — Wpd*(X,pQ). We rewrite the phase function 

x-P-e*(x,p) = x-p -e* {x, P ) +(P - P ) ■ [ (i - t)d PP e* (P Q + t[P - P }) dt [P - P ] . 

^ v ' JO 



=9{X,X) 



:!9 
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Figure 11. Comparision of the Schrodinger projection density p$(h) and the molecular 
dynamics density /?md for M — 200 (left plot) and M — 6400 (right plot). 
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Figure 12. Convergence rate of (8.13) for the observables gi and gi as defined in (8.14) 



The relation 



defines the function Y(P), and its inverse P(Y), so that the phase is a quadratic function in Y. The 
stationary phase expansion of an integral takes the form, see |10j . 



Y ■ dp P 0(Po)Y = (P-P ). (1 - t)d PP 9 (P + t[P - P }) dt [P - P Q ] 



w(P) e 



-iM 1/2 (X-P-8*(X,P)) 



dP 



(9.1) 



Vp9*(P )=X 



det 



d(P) 



d(X) 



1/2 



e ifsgn(3p P 0*(P o )) e -iM 1/2 e(X,X) 



00 A/T~k/'2 
k=0 \ l,j 



det 



d(P) 

W) 
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Algorithm 3 Algorithm for Example 2 



Input: An energy E, an one-dimensional potential function V, mass M, Schrodinger equation (8.7). 
Output: The Schrodinger projection density p$(h). 

1. Identify the right caustic point X + > satisfying X + = V^ 1 (E). For a fixed E € K, consider the 



periodic eigenvalue problem. Solve (8.7 1 numerically by constructing the discretised operator form of 



-(2M) x dxx + V using finite differences and denoted %W , and solve the eigenvalue problem 
3.15) H {h) P l = E.P, 

for the 10 eigenvalues closest to E using the Matlab eigenvalue solver eigs(H,10,E). Let Eq denote the 



eigenvalue closest to E and consider from now on solving (8.7) for the energy E and its corresponding 
eikonal equation \P 2 + V(X) = Eq. 

2. Determine 6*{P) by 



(P)= [ V P 9*(p)dp 
Jo 



3. Evaluate the Fourier integral 

3.16) u{X) := C * e iM^{-xp+e%p)) dp ^ ]x] > Xq ^ 

where Xq is chosen as the smallest value X > X + /2 such that is at a local maximum, and for 



\X\ < Xq compute u by (8.10) using the stationary phase method. 
4. Construct the approximate solution 

iCuiX^EQ-ViX))- 1 ^ \X\<X Q . 
\u(X)/y/\V x V(X)\ \X\>X Q . 

u{Xq){Eq-V{Xq)) 1 ^ 



$(X) := 

with 



C = 



^/\W x V(X )\u(Xq) 



5. 

Project $ onto the eigenspace {Ti}f =1 , J < 10 by Algorithm [2] to obtain a projection solution and 
compute its corresponding approximate density 

|$W|2(X) 
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